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1.  SUMMARY 


Regional  seismic  phases  (e.g.,  Pn,  Sn,  Pg,  Lg )  play  a  very  important  role  in  global 
monitoring  of  low  yield  underground  nuclear  tests.  Given  a  typical  paucity  of  teleseismic 
observations  for  a  small  event,  regional  phase  arrival  times  can  be  critical  for  locating  the 
source.  Numerous  empirical  observations  have  also  shown  that  regional  phases  hold  the 
keys  to  small  event  magnitude  and  yield  estimation  and  to  discrimination  between  small 
explosions  and  earthquakes  (e.g.,  Taylor  et  ah,  1989;  Kim  et  ah,  1997;  Walter  et  ah,  1995; 
Fisk  et  ah,  1996,  2005;  Taylor,  1996;  Taylor  and  Hartse,  1997;  Hartse  et  ah,  1997;  Patton, 
2001;  Xie,  2002;  Richards  and  Kim  2007;  Zhao  et  ah,  2008).  Use  of  travel  times  and 
amplitudes  of  regional  phases  for  location,  discrimination  and  yield  estimation  procedures 
requires  that  the  regional  velocity  and  attenuation  structures  be  sufficiently  well 
approximated  (e.g.,  Myers  et  ah,  2010).  As  a  result  there  have  been  extensive  efforts  to 
determine  regional  crustal  and  upper  mantle  velocity  models  and  tomographic  models  for 
attenuation  of  regional  phases.  These  efforts  commonly  proceed  independently,  with 
geometric  spreading  models  being  used  for  simple  or  generic  structures  even  while 
detailed  regional  velocity  models  are  determined  for  travel  times.  This  situation  is 
evolving,  with  an  increasing  emphasis  on  model-based  procedures  rather  than  on 
empirical  data  trends  entering  into  operational  systems,  and  this  presents  both  challenges 
and  opportunities  for  developing  self-consistent  model-based  procedures.  The  guiding 
theme  of  this  project  is  working  toward  unified,  self-consistent  correction  of  regional 
phases  for  amplitude  and  travel  time  effects  of  heterogeneous  Earth  structure.  Our  focus  is 
on  the  critical  regional  phase  Pn,  which  traverses  the  uppermost  mantle,  and  commonly  is 
assumed  to  have  very  simple  wave  behavior  controlling  its  amplitudes  and  travel  times.  In 
fact,  Pn  has  complex  sensitivity  to  Earth  structure  that  is  not  well  represented  by  standard 
processing  assumptions  (particularly  when  determining  attenuation  models  needed  for  the 
model-based  procedures),  motivating  use  of  regionally  appropriate  velocity  structures,  just 
as  are  motivated  for  precise  event  location.  Striving  toward  self-consistency  of  regional 
structures  used  for  travel  time  and  amplitude  analysis  is  a  logical  step  as  model-based 
approaches  are  adopted  in  the  operational  environment. 

2.  INTRODUCTION 

This  project  builds  on  a  prior  collaboration  in  which  we  developed  geometric 
spreading  corrections  for  regional  phases  Pn  and  Sn  for  simple  reference  velocity  models 
and  applied  these  to  develop  frequency-dependent  Pn  attenuation  models  for  Eurasia.  We 
confirmed  that  these  seismic  phases  have  acute  sensitivity  to  structure,  with  frequency- 
dependent  geometric  spreading  behavior  for  even  the  simplest  spherical  structures.  The 
sensitivity  to  mantle-lid  velocity  gradient  is  non-linear  and  pronounced,  but  can  be 
characterized  well  by  numerical  modeling  for  specific  structures.  Current  procedures  tend 
to  ignore  the  complex  geometric  spreading  of  these  phases  and  thus  likely  project  artificial 
frequency  dependence  into  attenuation  models  for  the  mantle  lid.  Practical  determination 
of  precise  regional  mantle-lid  velocity  gradients  is  challenging,  but  possible  if  curvature  of 
Pn  and  Sn  travel  time  branches  is  constrained  in  a  local  region.  This  is  exactly  what  is 
being  performed  in  RSTT  (regional  seismic  travel  time)  approaches,  although  data  scatter 
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and  trade-off  with  variable  crustal  legs  certainly  reduces  confidence  in  inferred  mantle-lid 
radial  velocity  gradients.  We  examine  the  predictions  of  RSTT  and  similar  models  for 
amplitude  effects,  considering  both  the  implications  for  associated  frequency-dependent 
attenuation  models  and  the  prospects  of  jointly  using  amplitude  and  travel  time 
information  to  better  constrain  the  regional  velocity  and  attenuation  structures.  Eurasian 
datasets  of  regional  Pn  phases  are  utilized,  and  our  previous  work  on  determining 
frequency-dependent  geometric  spreading  as  a  function  of  ID  velocity  model  variations  is 
expanded  to  provide  corrections  for  laterally  varying  regional  composites.  We  extending 
our  current  (up  to  10  Hz)  2D  finite-difference  modeling  efforts  to  more  realistic  media. 
The  ultimate  goal  is  to  work  toward  procedures  for  determining  regional  attenuation  and 
velocity  models  that  are  useful  for  self-consistently  correcting  both  travel-time  and 
amplitude  measurements  for  improved  event  location,  event  identification  and 
magnitude/yield  determination. 


3.  BACKGROUND  AND  TECHNICAL  APPROACH 

3.1  Pn  and  Sn  geometric  spreading 

Accurately  accounting  for  geometric  spreading  is  critical  for  the  development  of 
meaningful  regional-phase  attenuation  models.  This  is  particularly  true  for  Pn  and  Sn 
waves  because  the  nature  of  their  wave  propagation  renders  them  acutely  sensitive  to 
uppermost  mantle  velocity  structure  and  Earth’s  sphericity.  Even  simple  one-dimensional 
(ID)  velocity  models  can  produce  geometric  spreading  of  Pn  and  Sn  that  is  strongly 
dependent  on  frequency  (e.g.,  Sereno  and  Given,  1990;  Yang  et  ah,  2007;  Avants  et  al., 
2011).  If  frequency  dependence  of  the  geometric  spreading  is  neglected,  any  inferred 
attenuation  model  based  on  using  that  geometric  spreading  will  acquire  incorrect 
frequency  dependence. 

In  our  previous  modeling  efforts,  we  generated  synthetic  seismograms  for  a  ID, 
spherical  Earth  model  (Yang  et  al.,  2007;  Avants  et  al.,  2011)  to  simulate  Pn  and  Sn 
geometric  spreading  in  such  a  structure.  We  used  the  same  generic  spherical  Earth  model 
considered  by  Sereno  and  Given  (1990)  as  the  Base  Earth  Model  for  our  simulations 
(Figure  1)  and  developed  new  Pn  and  Sn  geometric-spreading  models  using  the  synthetics 
from  the  simulation. 
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Figure  1 .  Base  Earth  Model  used  for  Pn  and  Sn  simulations  and  the  development  of  new 
Pn  and  Sn  geometric-spreading  models.  Quality  factor  Q  is  infinite  throughout  the  model. 
Earth  flattening  transform  yields  the  ID  models  on  the  right,  which  have  mild  positive 
velocity  gradients  even  though  the  mantle  lid  has  constant  velocity  in  the  spherical 
structure. 

We  computed  synthetic  seismograms  using  the  source  parameterizations  and 
procedures  described  in  detail  by  Yang  et  al.  (2007)  and  extracted  the  Pn  and  Sn  portions 
of  the  synthetic  seismograms  using  fixed-velocity  windows.  We  measured  Pn  and  Sn 
amplitudes  from  the  spectra  calculated  using  the  Pn  and  Sn  synthetic  seismograms.  The 
10-Hz  Pn  amplitude  decay  in  the  spherical  Base  Earth  Model  is  shown  in  Figure  2.  Also 
plotted  in  the  figure  is  the  amplitude  decay  of  a  conical  head  wave  in  a  plane  one-layer- 
over-half-space  model  (Aki  and  Richards,  2002;  Eq.  6.26)  and  the  amplitude  decay  of 
infinite-frequency  direct  wave  in  a  spherical  Earth  model  from  ray  tracing.  At  distances 
close  to  the  critical  distance,  Pn  geometric  spreading  behaves  like  that  of  a  conical  head 
wave.  As  distance  increases,  Pn  spreading  starts  to  deviate  from  that  of  the  head  wave  and 
at  about  5°,  Pn  amplitudes  begin  to  increase.  In  the  range  between  the  critical  distance  and 
about  10°,  Pn  evolves  from  a  wave  similar  to  a  conical  head  wave  to  the  interference  head 
wave,  which  is  a  superposition  of  multiple  waves  reflected  from  the  Moho.  At  teleseismic 
distances,  the  direct-wave  spreading  approaches  that  of  the  infinite-frequency  diving  P 
wave  from  ray  tracing. 

Pn  geometric  spreading  in  a  spherical  Earth  model  is  not  only  different  from  that  of  a 
head  wave  as  is  shown  in  Figure  2,  but  also  frequency  dependent.  Figure  3  shows  the  Pn 
amplitude-variation  surface  as  a  function  of  distance  and  frequency  for  the  simple  Base 
Earth  Model.  The  strong  frequency  dependence  of  Pn  amplitudes  is  apparent.  Amplitudes 
at  higher  frequencies  are  affected  more  by  sphericity  than  are  lower-frequency  amplitudes. 
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Figure  2.  10-Hz  synthetic  Pn  amplitude  decay  in  the  spherical  Base  Earth  Model  with 
constant  mantle  velocities.  The  solid  line  depicts  the  theoretical  amplitude  decay  of  a 
conical  head  wave  in  a  plane  one-layer-over-half-space  Earth  model.  The  dashed  line  is 
the  amplitude  decay  of  infinite-frequency  direct  wave  in  a  spherical  homogeneous  Earth 
model  from  ray-tracing  calculations.  From  Yang  et  al.  (2007). 
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Figure  3.  Synthetic  Pn  amplitudes  as  a  function  of  epicentral  distance  and  frequency. 
From  Yang  et  al.  (2007). 
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Yang  et  al.  (2007)  presented  a  parameterized  description  of  Pn  geometric  spreading 
for  this  simple  reference  model.  The  amplitude  spectrum  of  Pn  can  be  parameterized  as 


A(r,ff,f)=K(f)M^0)G(r,f)exd  - 


(- 
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0t/> 


S(f) 


(1) 


with  the  new  geometric-spreading  model  expressed  as 
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In  Equation  1,  if  is  a  frequency-dependent  scaling  factor;  Mo  is  source  moment;  R  is 
source  radiation  pattern;  Q  is  Pn  quality  factor;  v  is  Pn  velocity;  S  is  receiver  site 
response;  r  is  epicentral  distance;  6  is  azimuth  angle  and  /  is  frequency,  ro  and  fo  are 
included  in  Equations  (2)  and  (3)  in  order  for  the  new  model  to  have  the  same  dimension 
as  standard  power-law  models  (e.g.,  Street  et  al.,  1975;  Sereno  et  al.,  1988).  Yang  et  al. 
(2007)  derived  the  coefficients  n,y  by  fitting  Equations  (2)  and  (3)  to  the  synthetic  data.  As 
a  result,  the  new  spreading  models  better  represent  Pn  and  Sn  geometric  spreading  in  a 
more  realistic  Earth.  The  main  differences  between  the  new  geometric-spreading  model 
and  the  standard  frequency-independent  power-law  model  are  the  addition  of  the  first  term 
in  the  exponent  and  the  frequency  dependence  of  parameters 

3.2  Pn  Spreading  Sensitivity  to  Mantle  Velocity  Gradients 

The  velocity  structure  of  the  mantle  lid  is  seldom  well  constrained.  The  most  direct 
constraint  is  provided  by  the  shape  of  the  Pn  travel  time  curve  as  a  function  of  distance 
over  a  wide-enough  distance  range  for  reliable  measurement,  and  even  when  there  is 
apparent  travel  time  branch  curvature,  it  is  difficult  to  distinguish  between  a  smooth 
gradient  with  depth  versus  constant  velocity  layering  with  small  step  increases.  Receiver 
function  methods  and  surface  wave  dispersion  inversions  sometimes  indicate  gradients  in 
the  lid  structure,  but  resolution  tends  to  be  poor,  especially  for  R-velocity.  As  a  result, 
reference  velocity  structures  often  assume  constant  or  near-constant  velocity  in  the  mantle 
lid  by  default,  so  the  Base  Earth  Model  is  not  a  bad  starting  point  in  many  cases. 
However,  there  are  regions  where  the  lid  velocity  appears  to  increase  or  decrease  with 
depth.  The  effects  of  positive  /-’-velocity  gradients  in  the  mantle  lid  on  the  frequency- 
dependent  spreading  of  Pn  are  illustrated  in  Figure  4.  The  basic  shape  of  the  spreading 
behavior  is  preserved,  but  the  distance-behavior  shifts  systematically  to  shorter  distances, 
with  overall  higher  Pn  amplitudes  as  the  gradient  increases.  It  is  possible  to  determine 
appropriate  coefficients  for  a  spreading  representation  like  Equation  (2)  for  any  specific 
gradient,  but  this  is  only  warranted  if  there  are  a  priori  constraints  on  the  specific  model  in 
the  region  for  which  the  spreading  is  to  be  utilized. 
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Figure  4.  Pn  amplitude  versus  distance  curves  similar  to  Figure  3,  but  for  the  Base  Earth 
Model  (red  curves,  for  no  physical  gradient  in  the  lid)  and  for  two  models  with  mild 
positive  gradients  in  the  mantle  lid.  The  suite  of  curves  for  each  case  corresponds  to 
frequencies  ranging  from  0.75  Hz  (lowest  curve)  to  12  Hz  (highest  curve).  From  Avants  et 
al.  (2011). 

3.3  Pn  Spreading  Sensitivity  to  Lateral  Mantle  Velocity  Volumetric  Heterogeneity 

In  Avants  et  al.  (2011),  we  explored  the  effects  of  lateral  heterogeneity  in  the  mantle 
lid  velocity  structure  on  Pn  spreading  using  a  2D  fourth-order  finite-difference  code  (Xie 
and  Lay,  1994).  We  included  random  lateral  volumetric  velocity  fluctuations  characterized 
by  RMS  velocity  and  varying  horizontal  and  vertical  averaging  functions  with  different 
length  scales  in  the  Model.  This  allowed  us  to  consider  2D  random  heterogeneity  models 
like  those  shown  in  Figure  5. 
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Figure  5.  Visualizations  of  2D  velocity  model  with  (top)  0.5%  (rms)  ^-velocity 
heterogeneity  in  the  upper  one  hundred  kilometers  of  the  mantle  lid,  and  (bottom) 
exponential  3%  (rms)  roughness  of  the  Moho  with  40  km  scale  averaging  length.  The 
background  model  is  the  EFT  version  of  the  Base  Earth  Model,  so  there  is  a  slight  positive 
velocity  gradient  across  the  crust  and  mantle.  From  Avants  et  al.  (2011). 


The  2D  finite-different  modeling  approach  used  at  that  time  could  not  achieve  the  very 
high  frequencies  of  the  ID  wave-number  integration  method,  so  the  synthetics  are  limited 
to  about  1  sec  dominant  period.  Analysis  of  effects  of  random  heterogeneity  requires  a 
statistical  sampling  of  the  effects  associated  with  different  realizations  of  the  random 
velocity  parameters.  Many  models  have  been  run  for  various  statistical  properties  of  the 
mantle  lid  structure,  with  complete  synthetic  seismograms  being  computed. 

Figure  6  shows  the  synthetic  ~1  Hz  Pn  amplitudes  for  ensemble  averages  of  5 
realizations  each  for  variable  heterogeneity  aspect  ratios  and  rms  velocity  fluctuations  for 
structures  like  that  at  the  top  in  Figure  5,  directly  compared  to  the  result  for  the  Base  Earth 
Model. 
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Figure  6.  Ensemble-averaged  ~1-Hz  Pn  amplitudes  as  a  function  of  distance  for  all 
configurations  of  random  volumetric  mantle  heterogeneity  simulated.  Pn  amplitude  for  the 
BEM  (diamonds)  are  compared  to  models  with  0.5%  RMS  (squares),  1.0%  RMS 
(triangles),  and  2%  RMS  (X's)  Vp  fluctuation  for  horizontal  (Ax)  and  vertical  (Ay) 
averaging  lengths  of  (a)  Ax  =  10  km,  Ay  =10  km  (isotropic);  (b)  Ax  =  20  km,  Ay  =  10 
km;  (c)  Ax  =  20  km,  Ay  =  6  km;  and  (d)  Ax  =  40  km,  Ay  =  3  km.  Each  value  is  the 
average  of  5  realizations  with  different  seed  kernels. 

An  interesting  behavior  in  Figure  6  is  that  the  presence  of  velocity  fluctuations  affects 
the  basic  shape  of  the  spreading,  and  in  some  cases  it  could  be  represented  well  by  a 
power-law  type  behavior  at  this  frequency  (~  1  Hz).  This  suggests  that  the  ID  spreading 
behavior  is  rather  delicate,  with  heterogeneity  disrupting  the  specific  interference  that 
gives  rise  to  the  complex  shape.  As  the  strength  of  the  velocity  functions  increases,  there 
is  progressive  reduction  of  curvature  of  the  amplitude-distance  trend,  and  the  overall 
amplitude  behavior  becomes  increasingly  power-law  like.  We  emphasize  that  these 
simulations  are  2D,  are  very  band-limited  so  any  frequency  dependencies  are  not  defined, 
and  have  prescribed  velocity  heterogeneity  spectra.  Whether  they  represent,  or  even 
approach  the  prevailing  situation  in  the  mantle  lid  is  an  open  question  and  constraints  can 
only  come  from  observed  data.  In  addition,  more  accurate  modeling  will  require 
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broadband  2D  and  3D  simulations.  Nevertheless,  the  results  demonstrate  the  effects  of 
strong  velocity  heterogeneity  on  Pn  geometric  spreading. 

3.4  Pn  Spreading  Sensitivity  to  Lateral  Moho  Topography 

While  volumetric  heterogeneity  is  expected  to  result  from  geological  processes  during 
crustal  formation  and  evolution,  the  same  is  true  for  Moho  irregularities  on  multiple 
scales.  Pn  and  Sn  interactions  with  the  Moho  near  the  source  and  receiver  as  well  as 
throughout  the  propagation  and  whispering  gallery  development  provide  sensitivity  to  the 
Moho  roughness.  Avants  et  al.  (2011)  parameterized  2D  exponential  (Figure  5)  and 
Gaussian  statistical  irregularities  in  the  Moho  depth  to  explore  the  effects  of  random 
structures,  along  with  simple  step-like  structure  in  the  Moho. 

The  exponential  model  is  much  richer  in  small-scale  structure  and  produces  longer 
enduring  coda  than  the  smoother  Gaussian  model,  but  the  overall  behavior  of  the  models 
is  similar  for  the  1  sec  dominant  period  (Figure  7).  For  an  80  km  horizontal  averaging 
function  the  amplitude  decay  effects  become  significant  for  5%  RMS  heterogeneity  in  the 
depth  of  the  Moho.  There  is  only  minor  sensitivity  to  the  choice  of  horizontal  averaging 
function  for  values  larger  than  40  km  for  a  given  statistical  model. 
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Figure  7.  Ensemble-averaged  ~1 -Hz  Pn  amplitudes  plotted  as  a  function  of  distance  for  all 
configurations  of  random  Moho  topography  heterogeneity,  grouped  by  percent  RMS  depth 
fluctuation.  Pn  amplitude  for  the  BEM  (diamonds)  are  compared  to  horizontal  averaging 
length  scales  of  20  km  (squares),  40  km  (triangles),  80  km  (X's),  and  160  km  (asterisks), 
for  percent  RMS  depth  fluctuations  of  (a)  1%;  (b)  3%;  and  (c)  5%.  From  Avants  et  al. 
(2011). 

3.5  Pn  Spreading  Model  Constrained  with  Observed  Data 

Yang  (2011)  evaluated  the  performance  of  Equation  (2)  in  correcting  observed  Pn 
amplitudes  in  Asia.  Figure  8  shows  the  1-Hz  Pn  amplitudes  after  source  and  geometric- 
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spreading  correction  using  Equation  (2).  The  figure  shows  that  Equation  (2)  provides  good 
spreading  correction  (leaving  a  linear  trend  in  log  amplitude  versus  distance)  only  over  a 
limited  distance  range  between  about  860  and  1400  km.  Outside  this  range,  corrected 
amplitudes  exhibit  undesirable  variations  (flattening  and  roughness)  relative  to  a  constant 
decay  rate  that  could  be  expected  for  average  attenuation  of  Pn. 
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Figure  8.  1-Hz  Pn  amplitudes  after  source  and  geometric-spreading  correction  using  the 
model  developed  by  Yang  et  al.  (2007).  The  red  line  is  the  400-point  moving  average 
showing  the  decay  trend  of  the  data. 

To  better  represent  Pn  spreading,  Yang  (2011)  developed  a  method  to  construct 
observation-based  and  region-specific  Pn  spreading  models.  In  his  method,  the  average  Pn 
Q  at  different  frequencies  is  calculated  from  source  and  spreading  (using  Equation  (2)) 
corrected  Pn  amplitudes.  These  Q  values  are  then  used  to  correct  raw  Pn  amplitudes  for 
attenuation.  The  resulting  amplitude  residuals  are  then  fit  with  an  observation-constrained 
spreading  model.  The  Pn  spreading  model  that  Yang  (2011)  developed  for  Asia  performs 
much  better,  resulting  in  a  set  of  spreading-corrected  amplitudes  that  have  a  constant 
distance  decay  slope  for  a  wide  distance  range  from  150  km  to  1600  km  (Figure  9). 
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Figure  9.  1-Hz  Pn  amplitudes  after  source  and  geometric-spreading  correction  using  the 
new  model  developed  by  Yang  (2011).  The  red  line  is  the  400-point  moving  average 
showing  the  decay  trend  of  the  data.  The  yellow  line  is  the  same  as  the  red  line  in  Figure 
8. 

3.6  The  Effect  of  Geometric  Spreading  on  Attenuation  Determination 

The  tradeoff  between  geometric-spreading  and  attenuation  makes  it  difficult  to 
accurately  estimate  one  parameter  if  the  other  parameter  is  not  well  constrained.  The 
problem  is  particularly  acute  for  Pn  and  Sn  due  to  reasons  discussed  above.  As  a  result,  an 
accurate  spreading  relationship  must  be  estimated,  even  for  individual  paths,  in  order  to 
develop  reliable  2D  attenuation  models  used  in  amplitude  corrections  for  discrimination 
and  yield  estimation. 

To  demonstrate  the  advantage  and  utility  of  a  physics-based  spreading  model  over 
traditional  power-law  model,  Figure  10  compares  the  measured  1-Hz  Pn  amplitudes  with 
the  spreading  model  of  Yang  et  al.  (2007)  and  a  power-law  model.  The  figure  shows  that 
the  Yang  et  al.  (2007)  model  has  a  consistently  slower  decay  rate  than  the  data  toward 
long  distances,  which  is  desirable  if  progressive  attenuation  occurs  to  produce  the  data 
trend  over  a  large  range  of  distance  (yellow  line).  The  power-law  model,  on  the  other 
hand,  decays  faster  than  the  data  toward  long  distances,  which  would  result  in  a  negative 
Q  when  trying  to  match  the  data  behavior. 

Figure  1 1  demonstrates  this,  showing  examples  of  Q  estimates  obtained  using  either 
the  model  of  Yang  et  al.  (2007)  or  the  power-law  model  for  multiple  frequencies.  The 
model  of  Yang  et  al.  (2007)  yields  reasonable  Q  values  at  all  frequencies  whereas  using 
the  power-law  model  results  in  either  negative  or  extremely  large  Q  values.  These 
examples  show  that  even  though  the  model  of  Yang  et  al.  (2007)  is  not  optimal,  it 


11 

Approved  for  public  release;  distribution  is  unlimited. 


provides  better  spreading  corrections,  due  to  its  physics  basis,  than  a  power-law  model, 
which  corresponds  to  no  known  Earth  structure. 

—2 - 1 - 1 - 1 - 1 - 1  i  i  i  i 


-3 


-4 

-5 


-7 

-8 

-9 

-10 

-11 

-12 

0  200  400  600  800  1000  1200  1400  1600 

Figure  10.  Comparison  of  1-Hz  Pn  amplitude  decay  with  decays  predicted  by  a  power-law 
spreading  model  and  by  the  new  spreading  model  of  Yang  et  al.  (2007).  The  yellow  line 
is  the  400-point  moving  average  of  the  Pn  data  points  to  more-clearly  delineate 
the  amplitude-decay  trend. 

The  foregoing  review  of  accomplishments  and  enhanced  understanding  of  geometric 
spreading  for  Pn  (all  attributes  hold  for  Sn  as  well)  indicates  that  reliable  models  of 
attenuation  can  be  obtained  only  if  reliable  geometric  spreading  is  computed  for  specific 
regional  structure.  Thus,  a  strategy  for  determining  regional  structures  with 
observationally-constrained  mantle-lid  velocity  gradients  is  needed. 

The  strategy  of  regional  seismic  travel  time  (RSTT)  model  parameterization  advanced 
by  Myers  et  al.  (2010)  adopts  a  geographic  node  description  of  crustal  and  upper  mantle 
models,  where  the  mantle  portion  is  parameterized  by  a  linear  P-wave  velocity  gradient. 
This  parameterization  is  very  similar  to  that  of  the  Base  Earth  Model  used  in  the  Pn  and  Sn 
geometric  spreading  calculations  of  Yang  et  al.  (2007)  and  Avants  et  al.  (2011)  (Figure 
12).  Myers  et  al.  (2010)  used  regional  Pn  travel  times  to  develop  a  tomographic  model  of 
Pn  velocity  variation  and  mantle-gradient  variation  in  the  mantle  lid  (Figure  13).  While  the 
parameterization  is  not  explicitly  intended  to  predict  amplitude  behavior,  nor  are  the 
resulting  velocity  models  believed  to  be  accurate  in  each  case,  the  models  can  be  used  as  a 
starting  point  to  explore  amplitude  behavior. 


1-Hz  Pn  amplitude 
400-point  moving  average 
power-law  spreading 


new  spreading 
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new  spreading  model 
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Figure  1 1 .  Comparison  of  average  Q  estimates  from  Pn  amplitudes  corrected  for  source 
and  geometric-spreading  effects  using  different  spreading  models.  Q  estimates  at  1,  4  and 
8  Hz  are  compared. 


(a) 


Geographic  node  sampling 


(b)  Model  parameterization 
at  each  node 


Velocity  (km/s) 

5  6  7  8  9 


Figure  12.  Global  model  parameterization,  (a)  An  example  tessellation  with  approximately 
5°  grid  spacing.  The  inset  shows  the  1°  used  in  the  study  by  Myers  et  al.  (2010).  Color 
indicates  Moho  depth  of  the  starting  model,  (b)  An  example  velocity  vs.  depth  profile  as 
defined  at  each  node.  The  mantle  portion  of  the  profile  is  specified  by  the  velocity  at  the 
crust/mantle  interface  and  a  linear  gradient.  From  Myers  et  al.  (2010). 


Regional  phase  arrival  times  were  used  by  Myers  et  al.  (2010)  to  develop  a 
tomographic  model  (Figure  13)  for  variation  in  velocity  of  the  uppermost  mantle  and 
variation  in  the  strength  of  the  mantle  gradient.  This  provides  regional  models  that  fit 
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travel  times,  but  the  precision  of  the  velocity  gradients  in  the  mantle  is  impacted  by  scatter 
in  the  data,  averaging  of  nodes  in  a  given  region,  and  lack  of  resolution  of  curvature  of  Pn 
travel  times.  While  the  parameterization  is  not  explicitly  intended  to  predict  amplitude 
behavior,  the  models  can  be  used  to  do  so. 

We  examine  RSTT  derived  models  and  associated  travel  time  data,  with  corresponding 
regional  phase  amplitude  measures  as  a  function  of  frequency  to  evaluate  use  of  travel 
time  models  as  first-order  structures  for  computing  regional  geometric  spreading 
corrections  (thereby  having  unified  model  frameworks  for  travel  time  and  amplitudes),  and 
as  a  starting  point  for  revision  of  the  models  guided  by  observed  amplitude  behavior  that  is 
not  well  matched. 
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Figure  13.  Comparison  of  starting  and  RSTT  models,  (a)  Velocity  below  the  Moho  for 
starting  model  and  (b)  RSTT  model,  (c)  Mantle  gradient  (km/sec/km)  for  starting  model 
and  (d)  RSTT  model.  From  Myers  et  al.  (2010). 
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4.  RESULTS 


Building  on  the  back-ground  results  discussed  above  from  an  earlier  contract,  we  have 
performed  additional  forward  modeling  of  Pn  geometric  spreading  in  various  ID  and  2D 
structures,  analyzed  regional  variations  in  Pn  around  the  North  Korean  test  site  and 
modeled  effects  of  oceanic  path  segments,  and  explored  the  RSTT  models  and  empirical 
spreading  correction  approach  for  Eurasia,  as  discussed  in  this  section. 

4.1  One-Dimensional  Velocity  Models 

To  simulate  Pn-wave  propagation,  a  set  of  ID  and  2D  velocity  models  with  different 
velocity  gradients  were  designed.  These  models  are  illustrated  in  Figures  14  and  15.  All 
model  parameters  shown  are  functions  of  radius  before  earth  flattening  transformation 
(EFT).  A  constant- velocity  lid  model  (Const-Lid)  is  shown  in  Figure  14  and  listed  in 
Table  1.  For  comparison,  we  also  calculated  Pn  propagation  in  the  Const-Lid  model 
without  applying  the  EFT  (Const-Lid  w/o  EFT).  A  velocity  model  with  lid  radial  gradient 
of  lxlO”3  sec”1  (Gradient-0.001)  is  shown  in  Figure  14  and  described  in  Table  2.  Velocity 
models  with  0x10  3  sec”1  and  2xl0”3sec“1  gradients  (Gradient-0.000  and  Gradient-0.002) 
are  obtained  similarly.  These  models  connect  to  a  steep  IASP-91  gradient  at  greater  depth. 

Vp 

Ln  <Tl  >vj  00  tD 

Lno'iLn^LnooLntisLn 


Figure  14.  One-dimensional  velocity  models  used  for  Pn  computation.  The  structures  are 
spherical  models  prior  to  Earth  Flattening  Transform  (EFT). 
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Figure  15.  Velocity  models  for  Gradient-0.003  and  Gradient-0.004.  Note,  here  the  P- 
velocity  beneath  the  Moho  is  7.8  km/s,  so  that  the  strong  positive  gradients  do  not  lead  to 
un-physically  high  P  velocities. 


Table  1.  Constant  velocity  lid  mode 

,  Const-Lid  (radia 

parameters  before  EFT) 

Depth  of  the 
layer 

Vp  (km/s) 

Vs  (km/s) 

P  ( g  cm-3 ) 

crust 

40  km 

6.5 

3.75 

2.7 

mantle 

infinity 

8.1 

4.45 

3.32 

Table  2.  Velocity  model  with  1x10  3  sec  1  gradient,  Gradient-0.001  (radius  parameters 
before  EFT) _ 


Vp  (km/s) 

Vs  (km/s) 

p  (g/cmA3) 

crust 

40  km 

6.5 

3.75 

2.7 

mantle 

Starting  from  8. 1 
km/s  below 

Moho, 

increasing  with  a 
lxlO”3  sec”1 
gradient  until 
reach  to  a  depth 
where  Vp  has  the 
same  values  as 
model  IASP-91 
thereafter. 

Converted  from 
Vp  adopting 
Vp/Vs  values  in 
IASP-91 

From  PREM 
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4.2.  Laterally  Varying  Velocity  Models 


Velocity  models  with  laterally  varying  velocity  gradients  are  used  for  testing  Pn-wave 
propagation  in  these  environments.  These  models  are  composed  of  two  sections  with 
different  gradients.  One  section  from  0  -  450  km,  another  from  550  km  -  1000  km. 
Between  them,  from  450  km  -  550  km,  the  linear  interpolation  is  used.  These  velocity 
models  include: 

Gradient-0.000-0.001  (velocity  gradient  varies  from  0xl0”3sec_1  to  lxlO”3  sec”1), 

Gradient-0.00 1  -0.000  (velocity  gradient  varies  from  1  x  1 0”3  sec”1  to  0x1 0“3  sec”1 ), 

Gradient-0.000-0.002  (velocity  gradient  varies  from  0xl0”3sec_1  to  2xl0”3sec”1), 

Gradient-0.002-0.000  (velocity  gradient  varies  from  2  x  10”3  sec”1  to  0  x  10“3  sec”1). 

Figure  16  shows  displays  of  these  models  relative  to  model  Gradient-0.000. 

Distance  (km) 


□  250  SOD  75C  10OC 


O.ttKME+PO  <l.1?4fi£+06 


Figure  16.  Laterally  varying  velocity  models.  Shown  here  are  gradient  models  subtracting 
Gradient-0.000.  From  top  to  bottom  are  for  Gradient-0.000-0.001,  Gradient-0.00 1-0.000, 
Gradient-0.000-0.002,  and  Gradient-0.002-0.000,  respectively. 
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4.3  Pn  Simulations  for  ID  Velocity  Models 


Extending  the  bandwidth  of  finite-difference  calculations  for  large-distance  Pn 
synthetics  was  one  of  our  initial  priorities.  The  fourth-order  code  used  by  Avants  et  al. 
(2011)  was  adapted  to  and  tested  on  our  cluster,  and  stable  P-SV  and  SH  calculations  were 
achieved  for  ranges  of  up  to  1000  km  and  frequencies  up  to  10  Hz,  for  models  with  a  flat 
ffee-surface  and  laterally  varying  structure,  either  with  or  without  small-scale  statistical 
heterogeneity  being  included.  This  enabled  examination  of  Pn  geometric  spreading 
behavior  in  laterally  varying  2D  structure  over  broader  bandwidth  than  in  prior  work. 
Table  3  indicated  parameters  used  in  the  finite-difference  calculations. 


Table  3.  Parameters  used  for  FD  calculations 


Model  size  (horizontal) 

1000  km 

Model  size  (Vertical) 

400  km 

dx 

0.05  km 

dz 

0.05  km 

grid  size 

20003  x  8003 

dt 

3.5  ms 

elapsed  time 

145  sec. 

time  steps 

41488 

Source  depth 

15  km 

Type  of  source 

Explosion  source 

Figure  17  shows  wavefield  snapshots  for  the  ID  structure  Const-Lid  w/o  EFT  out  to 
ranges  of  900  km.  In  this  case  there  is  no  induced  positive  velocity  gradient  from  the 
earth-flattening,  so  the  model  is  the  classic  constant  velocity  layered  structure  that  gives 
rise  to  a  true  head  wave  with  Pn  amplitude  decay  following  the  solid  line  in  Figure  2.  The 
wavefield  behavior  is  quite  different  when  the  EFT  is  included,  as  seen  in  Figure  18.  The 
Pn  amplitudes  are  now  significantly  higher  at  large  range.  This  corresponds  to  the  basic 
behavior  described  in  the  previous  section,  with  turning  ray  behavior  of  the  multiples  in 
the  lid  modifying  the  geometric  spreading  of  the  first  arrival.  Increasing  the  positive 
gradients  beyond  the  mild  gradient  from  EFT  further  increases  the  Pn  amplitudes,  as 
shown  in  the  snapshots  for  model  Gradient-0.001  shown  in  Figure  19  and  for  model 
Gradient-0.002  in  Figure  20.  Figure  21  compares  the  waveforms  for  these  four  models  at 
a  time  lapse  of  120  s.  The  amplitude  normalization  is  the  same  for  all  four  models  so  that 
amplitudes  are  directly  comparable.  Note  the  strengthening  of  the  Pn  amplitudes 
proportional  to  velocity  gradient  and  the  clear  development  of  underside  reflections  from 
the  Moho  as  the  gradient  increases.  These  calculations  replicate  the  findings  from  our 
earlier  project,  but  extend  the  frequency  band  of  the  finite-difference  calculations  to  about 
10  Hz  for  distances  of  1000  km,  as  had  only  been  achieved  with  frequency- wavenumber 
integration  previously.  This  enables  exploration  of  high  frequency  Pn  geometric  spreading 
behavior  in  2D  models  for  the  first  time. 
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TIME  =  40.00  sec 
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distance  (km) 


700  800  900  1000 


Figure  17.  Snapshots  of  Pn  waves  propagating  in  the  Const-Lid  w/o  EFT  model  using  the 
enhanced  2D  finite-difference  code  implementation. 
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Figure  18.  Snapshots  of  Pn  waves  propagating  in  the  Const- Lid  model.  Note  the  enhanced 
Pn  amplitude  resulting  from  the  positive  velocity  gradients  from  the  EFT  (compare  to 
Figure  17). 
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Figure  19.  Snapshots  of  Pn  waves  propagating  in  velocity  model  Gradient-0.001  with 
1  x  10-3  sec-1  velocity  gradient  in  the  upper  mantle. 
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Figure  20.  Snapshots  of  Pn  waves  propagating  in  velocity  model  Gradient-0.002  with 
2  x  10“3  sec-1  velocity  gradient  in  the  upper  mantle. 
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Figure  21.  Comparison  between  snapshots  in  models  with  different  velocity  gradients  (t  = 
120  sec).  For  easy  comparison,  the  same  normalization  factors  are  used  for  all  models. 
From  top  to  bottom  are  snapshots  in  model  Const-Lid  w/o  EFT,  Const-Lid,  Gradient- 
0.001  and  Gradient-0.002,  respectively.  Note  the  different  Pn  amplitude  and  different 
mantle  P-wave  wavefront  expansion.  In  models  with  large  gradient,  the  wave  front 
bouncing  from  the  bottom  of  the  Moho  can  be  seen. 

Figure  22  shows  amplitudes  of  the  Pn  arrivals  for  the  same  four  models  in  reduced- 
velocity  (8.1  km/s)  distance  profiles.  The  amplitudes  for  the  two  Const-Lid  models  are 
enhanced  by  a  factor  of  3  relative  to  those  in  Gradient  models.  The  increase  in  Pn 
amplitude  with  increasing  range  and  earlier  arrival  Pn  arrival  for  the  Gradient-0.001  and 
Gradient-0.002  models  demonstrates  the  curvature  of  the  geometric  spreading  curve  seen 
in  Figure  4. 
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(a)  Const  lid  w/o  EFT 


(b)  Const  lid 


(c)  Gradient  0.001 


(d)  Gradient  0.002 


Figure  22.  Comparison  between  Pn  waves  calculated  in  various  ID  velocity  models,  with 
amplification  factors  used  in  (a)  and  (b)  being  3  times  those  used  for  (c)  and  (d).  As 
references,  the  solid  lines  are  group  velocities  of  8.3,  8.2,  8.1,  8.0  and  7.9  km/s.  The  center 
mark  indicates  the  picked  arrival  time  and  the  two  short  marks  indicate  the  sampling 
windows  for  Pn  amplitude  calculation.  However,  the  synthetic  data  are  quite  different 
from  the  real  data.  The  synthetics  are  more  impulsive,  while  the  real  data  are  more  like 
wave  trains.  A  reduced  time  t-r/8.1  is  used  for  horizontal  coordinate. 
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For  even  stronger  positive  gradients,  as  in  models  Gradient-0.003  and  Gradient-0.004, 
the  Pn  amplitude  and  waveform  continue  to  change  and  the  sub-Moho  multiples  become 
increasingly  apparent.  This  is  shown  by  the  waveform  profiles  in  Figure  23. 


(a)  Gradient  0.003 


(b)  Gradient  0.004 


Figure  23.  Comparison  between  Pn  waves  calculated  in  velocity  models  Gradient-0.003 
and  Gradient-0.004.  Note,  due  to  the  slow  upper  mantle  P-wave  velocity  (7.8  km/s  below 
the  Moho),  the  Pn  velocity  is  low.  For  reference,  the  solid  lines  are  group  velocities  of  8.0, 
7.9,  7.8,  7.7  and  7.6  km/s.  A  reduced  time  t-r/7.8  is  used  for  the  horizontal  coordinate. 

Figure  24  shows  Pn  arrival  time  picks  for  broadband  synthetics  for  4  of  the  models.  In 
the  Const-Lid  w/o  EFT,  the  Pn  arrival  time  has  almost  constant  move-out.  For  the  3 
models  with  positive  gradients,  including  the  contribution  from  the  EFT,  Pn  arrives  with 
higher  apparent  velocity  at  larger  distances  (also  refer  to  the  waveform  figures  and  group 
velocity  lines  in  Figure  22).  These  arrival  times  are  picked  using  a  very  simple  method 
and  no  frequency  filter  is  applied.  Slightly  different  arrival  time  picks  will  be  made  for 
different  passbands.  The  travel  times  also  include  a  time  shift  from  the  source  time 
function  in  the  synthetics.  This  figure  demonstrates  that  the  travel  time  variations  from 
varying  lid  velocity  gradients  can  be  observed,  and  if  a  ID  velocity  model  with 
lithospheric  positive  velocity  gradient  is  appropriate  for  a  given  region,  it  could  be 
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constrained  by  travel  times  and  used  to  develop  self-consistent  geometric  spreading  curves 
for  Pn  and  Sn. 
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Figure  24.  Arrival  time  picks  from  profiles  of  synthetics  in  Figure  22  displaying  the  weak 
curvature  introduced  by  linear  mantle  lid  gradients.  While  the  variations  are  subtle,  given 
sufficient  observations  and  laterally  uniform  structure,  the  mantle  gradient  can  be 
constrained  by  Pn  arrival  time  picks,  providing  velocity  models  that  can  be  used  for  self- 
consistent  regional  amplitude  spreading  calculations. 


4.4  Pn  Geometric  Spreading  for  ID  and  2D  Velocity  Models 

For  the  ID  models  in  Figures  14  and  15,  we  computed  frequency-dependent  geometric 
spreading  amplitude  decay  out  to  1000  km  using  the  finite-different  calculations, 
correcting  for  excitation  of  a  point-source  versus  a  line-source.  This  yields  the  amplitude- 
distance  curves  shown  in  Figure  25.  These  patterns  are  very  similar  to  wavenumber- 
integration  calculations  (Figure  4)  for  the  ID  velocity  models  computed  in  Avants  et  al. 
(2011).  The  complex  frequency-dependence  of  Pn,  with  minima  in  the  amplitude-distance 
curves  and  increasing  amplitudes  at  larger  distances,  particularly  for  higher  frequencies, 
being  evident  for  all  structures  with  a  positive  gradient  in  the  ID  flat  earth  structure. 
Comparison  of  the  amplitude  behavior  for  these  models  for  different  narrow  frequency 
bands  is  displayed  in  Figure  26.  For  lid  gradients  in  excess  of  0.002  sec'1,  the  distance  and 
frequency  dependence  are  quite  dramatic  as  shown  in  Figure  27  and  Figure  28,  with 
complex  plateauing  and  oscillation  of  the  high  frequency  components.  Some  of  the  large 
distance,  high  frequency  behavior  is  not  well-resolved,  but  very  similar  behavior  is  found 
for  the  completely  independent  frequency-wavenumber  integration  synthetics  in  Figure  4. 
Overall,  these  calculations  provide  a  good  level  of  confidence  in  the  broadband  Pn 
synthetics,  and  we  proceed  to  consider  results  for  2D  models. 
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Figure  25.  Frequency-dependence  of  Pn  amplitude  vs.  distance  relations  on  log- log  scale, 
(a)  to  (d)  are  for  different  velocity  models.  The  frequency  filters  are  indicated  by  the 
symbols.  Rms  amplitudes  are  used  in  this  figure,  and  all  curves  are  normalized  at  300  km. 
As  expected,  in  the  Const.-Lid  w/o  EFT,  the  amplitude-distance  relation  for  a  true  head 
wave  is  frequency  independent. 

The  Pn  amplitude  distance  curves  for  2D  models,  with  lateral  variations  in  the  radial 
gradients,  as  displayed  in  Figure  16  are  shown  in  Figures  29-32.  Figure  29  shows  the 
variation  between  end-member  ID  models  Gradient-0.000  and  Gradient-0.001,  with  the 
near  source  environment  and  receiver  environment  interchanging  between  these  two 
models.  All  models  have  positive  gradients  from  the  EFT.  Note  the  difference  in 
geometric  spreading  for  sources  in  different  near-source  velocity  structure.  This  indicates 
that  it  is  not  straightforward  to  use  a  single  average  model  to  represent  the  geometric 
spreading  in  a  laterally  varying  model.  The  lateral  averaging  of  the  2D  structures  is  such 
that  the  closer  distances  have  amplitude  variations  consistent  with  the  ID  end-member 
structure  near  the  source,  while  the  larger  distances  have  geometric  spreading  similar  to 
the  ID  end-member  structure  near  the  receiver.  This  is  emphasized  in  Figure  30,  for  the 
2. 5-5.0  Hz  passband,  where  the  ID  and  2D  spreading  curves  are  superimposed.  The  2D 
spreading  behavior  is  overall  intermediate  to  the  end-members,  but  does  depend  on 
direction  of  the  paths. 
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Figure  26.  Amplitude  vs.  distance  on  log-log  scale  for  different  models,  for  narrow-band 
filtered  signals  with  central  frequencies  and  passbands  of:  (a)  0.75  Hz  (0. 5-1.0  Hz),  (b) 
3.75  Hz  (2. 5-5.0  Hz),  and  (c)  7.0  Hz  (6.0-8.0  Hz). 
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(a)  Gradient  0.003 
3.0  T — 


Figure  27.  Frequency-dependence  of  amplitude  vs.  distance  relations  in  log-log  scale,  (a) 
In  Gradient-0.003  model,  and  (b)  Gradient-0.004  model.  Refer  to  Figure  15  for  velocity 
models. 


—  Gradient  0.003 
Gradient  0.004 


Figure  28.  Comparison  of  amplitude  vs.  distance  relations  for  Gradient-0.003  and 
Gradient-0.004  models  at  frequency  3.7  Hz. 
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Figure  29.  Amplitude-distance  relations  in  different  ID  and  2D  laterally  varying  velocity 
models  and  frequency  bands.  The  models  are:  (a)  Gradient-0.000,  (b)  Gradient-0.000- 
0.001,  (c)  Gradient  model  0.001-0.000,  and  (d)  Gradient-0.001.  Note  that  the  amplitudes 
for  6-8  Hz  and  8-10  Hz  bands  are  very  close. 
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Figure  30.  Comparisons  for  amplitude-distance  relations  for  ID  and  2D  laterally  varying 
models.  The  used  frequency  band  is  2. 5-5.0  Hz.  Note  the  transition  zone  in  the  2D 
structures  is  between  450  -  550  km. 
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Similar  comparisons  are  shown  for  the  models  involving  stronger  radial  gradients  in 
the  lid  in  Figures  31  and  32,  ranging  from  end-member  models  Gradient-0.000  to 
Gradient-0.002.  The  lateral  averaging  behavior  is  similar  to  the  Gradient-0.001  models, 
although  the  stronger  amplitude  variations  associated  with  the  strong  gradient  in  Gradient- 
0.002  enhance  the  variations  in  the  2D  models.  The  frequency-dependence  at  large  ranges 
is  strong  for  all  models,  but  the  distance  variation  in  the  amplitudes  is  distinct  for  each 
model,  as  emphasized  in  Figure  30,  which  compares  the  amplitude-distance  curves  for  the 
2. 5-5.0  Hz  band.  The  complexity  of  the  frequency-dependent  amplitudes  for  the  laterally 
varying  models  and  the  strong  direction-dependence  of  the  amplitude  behavior  would  be 
mapped  into  complex  frequency-dependence  of  attenuation  and  possible  anisotropic 
attenuation  if  the  effects  of  the  elastic  velocity  structure  are  not  correctly  accounted  for. 
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Figure  31.  Amplitude-distance  relations  in  different  laterally  varying  velocity  models  and 
frequency  bands.  The  models  are:  (a)  Gradient-0.000,  (b)  Gradient-0.000-0.002,  (c) 
Gradient  model  0.002-0.000,  and  (d)  Gradient-0.002.  Note  that  the  amplitudes  for  6-8  Hz 
and  8-10  Hz  bands  are  very  close. 
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Gradient  0.000 
Gradient  0.000  -  0.002 
Gradient  0.002  -  0.000 
Gradient  0.002 


Figure  32.  Comparisons  of  amplitude-distance  relations  for  different  laterally-varying 
velocity  models.  The  frequency  band  is  2. 5-5.0  Hz.  Note  the  transition  zone  in  the  2D 
models  is  between  450  -  550  km. 


4.5  Pn  Spreading  in  Models  With  Fine-Scale  Heterogeneity 


It  is  now  well  established  by  the  work  described  above  that  Pn  (and  Sn)  can  not 
usually  be  treated  as  a  simple  conical  headwave  developed  along  a  flat  Moho 
discontinuity.  Figure  33  schematically  depicts  the  factors  that  affect  the  amplitude  of  Pn. 
Starting  from  a  conventional  headwave  case,  where  a  constant  low-velocity  crust  (the 
upper  half-space)  is  welded  to  a  constant  high-velocity  mantle  (the  lower  half-space)  at  a 
planar  Moho  discontinuity  (Figure  33a),  a  source  located  in  the  crust  can  excite  waves 
propagating  from  the  source  to  the  refracting  point  at  the  Moho,  and  then  propagating 
horizontally  in  the  mantle  immediately  below  the  Moho.  The  energy  of  this  high-speed 
mantle  wave  can  leak  back  to  the  crust  and  be  observed  as  the  headwave  in  the  upper  half 
space.  The  upper  mantle  leg  is  composed  of  waves  with  imaginary  vertical  wavenumber. 
In  this  circumstance,  there  is  no  vertically  propagating  mode.  The  energy  immediately 
below  the  Moho  is  quickly  depleted  without  being  recharged  from  the  deeper  depth, 
resulting  in  rapid  decay  of  the  observed  headwave  amplitude.  Actual  Pn  wave 
propagation,  particularly  for  Pn  amplitude  variations  in  a  spherical  Earth,  is  much  more 
complicated.  Cerveny  and  Ravindra  [1971]  and  Hill  [1973]  theoretically  investigated  the 
behavior  of  Pn  waves  in  a  spherical  Earth  model  and  regarded  this  phase  as  an  interference 
of  multiple-diving  waves  reflected  from  the  underside  of  the  Moho  discontinuity.  Sereno 
and  Given  [1990]  studied  Pn  waves  in  flat  and  spherical  earth  models  and  found  that 
Earth’s  sphericity  alone  causes  a  significant  complexity  in  Pn  geometrical  spreading  and 
this  phenomenon  is  frequency  dependent.  The  sphericity  in  the  Earth  is  equivalent  to  a 
positive  velocity  gradient  in  the  mantle  (Figure  33b),  which  causes  whispering  gallery 
phases  reflected  from  the  underside  of  the  Moho.  In  this  way  the  energy  at  the  great  depth 
can  feed  the  depleted  energy  below  the  Moho,  making  the  observed  Pn  wave  much 
stronger  than  the  headwave.  In  addition,  in  a  simple  two  constant  velocity  half-space 
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system,  there  is  no  characteristic  length  scale.  Waves  of  different  wavelengths  (or 
frequencies)  behave  equally.  Therefore,  the  geometrical  spreading  function  for  headwave 
is  frequency  independent.  However,  in  a  model  with  upper  mantle  velocity  gradient,  the 
gradient  serves  as  a  characteristic  scale,  and  the  Pn  wave  geometric  spreading  function 
becomes  frequency  dependent. 

If  lateral  velocity  variations  exist  in  the  upper  mantle,  scattering  may  affect  the  Pn 
wave  propagation  in  a  way  that  is  quite  different  from  its  effects  on  other  waves  (Figure 
33c).  For  most  body  and  surface  waves,  the  scattering  process  deflects  the  wave  energy 
away  from  its  main  propagation  direction,  causing  the  wave  to  lose  energy.  This  process  is 
known  as  scattering  attenuation.  However,  the  Pn  wave  is  confined  by  the  Moho 
discontinuity  and  the  upper  mantle  velocity  gradient.  In  a  laterally  smooth  waveguide, 
relatively  little  energy  escapes  from  the  main  propagation  channel  to  radiate  upward  to 
form  the  observed  Pn.  The  presence  of  small-scale  heterogeneities  in  the  upper  mantle  can 
cause  scattering  and  its  effect  on  Pn  wave  propagation  are  twofold.  First,  as  for  other 
waves,  the  scattering  can  cause  attenuation  along  the  upper  mantle  channel.  Second,  the 
scattering  makes  an  efficient  mechanism  to  increase  the  observed  Pn  amplitude  by 
deflecting  the  energy  away  from  the  main  propagation  direction  in  the  upper  mantle. 
Compared  to  otherwise  very  weak  Pn  energy,  the  scattered  energy  can  be  strong,  causing 
prominent  enhancement  of  the  Pn  amplitudes.  This  often  causes  somewhat  of  a  dilemma 
for  the  Pn  wave,  as  to  whether  the  scattering  should  be  treated  as  scattering  attenuation,  or 
as  geometrical  spreading,  but  both  are  elastic  effects  controlled  by  the  velocity  structure. 


(a)  Velocity 

1 


Moho 


Observed  Pn 


(b) 

s 

(c) 

1 

Figure  33.  Pn  waves  in  different  crust-upper  mantle  velocity  models,  with  (a)  simple 
constant  velocity  crust  and  upper  mantle,  and  a  planar  Moho,  giving  rise  to  a  true  conical 
head  wave;  (b)  vertical  velocity  gradient  in  the  upper  mantle,  giving  rise  to  the  frequency- 
dependent  spreading  examined  earlier,  and  (c)  laterally  varying  fine-scale  velocity 
heterogeneity  in  the  upper  mantle. 
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In  this  section  we  introduce  fine-scale  heterogeneities  to  our  smoothly  varying  velocity 
gradient  models  to  explore  the  effect  on  the  Pn  geometric  spreading  behavior.  The  goal  is 
to  extend  the  calculations  of  Avants  et  al.  (2011)  to  higher  frequency  and  to  thereby 
evaluate  frequency-dependent  effects  of  the  fine-scale  structure.  The  random  velocity 
models  are  generated  by  adding  velocity  and  density  perturbations  to  the  mantle  part  of 
ID  velocity  models  illustrated  in  Figure  14.  These  models  have  a  constant  velocity  crust. 
Below  the  Moho  discontinuity,  the  mantle  has  a  given  vertical  velocity  gradient  of  0.000, 
0.001  or  0.002.  Below  this  gradient  layer,  the  mantle  follows  the  Iasp-91  earth  model.  The 
random  velocity  perturbations  have  an  exponential  spectrum,  and  are  characterized  by  the 
horizontal  and  vertical  correlation  lengths,  and  the  rms  velocity  perturbation.  Both  P  and  S 
wave  velocities  have  the  same  percentage  perturbation,  and  the  density  perturbation  is  half 
of  the  velocity  perturbation.  The  model  parameters  are  listed  in  Table  4. 

Table  4.  2D  Velocity  Models  With  Fine  Scale  Heterogeneity 


No 

Background  model 

Type  of 
random 
function 

Correlation 

length 

RMS 

perturbation 

Radom  seeds 
calculated 

1 

Const  lid  (0.000) 

exponential 

20  km  x  6  km 

0.5% 

3, 5, 7, 8, 9 

2 

Const  lid  (0.000) 

exponential 

20  km  x  6  km 

1.0% 

3, 5, 7, 8, 9 

3 

Const  lid  (0.000) 

exponential 

20  km  x  6  km 

2.0% 

3, 5, 7, 8, 9 

4 

Const  lid+0.001  grad. 

exponential 

20  km  x  6  km 

0.5% 

3, 5, 7, 8, 9 

5 

Const  lid+0.001  grad. 

exponential 

20  km  x  6  km 

1.0% 

3, 5, 7, 8, 9 

6 

Const  lid+0.002  grad. 

exponential 

20  km  x  6  km 

1.0% 

3, 5, 7, 8, 9 

7 

Const  lid+0.001  grad. 

exponential 

40  kmx  10  km 

0.5% 

3, 5, 7, 8, 9 

8 

Const  lid+0.001  grad. 

exponential 

40  kmx  10  km 

1.0% 

3, 5, 7, 8, 9 

9 

Const  lid  +  0.001  grad. 

exponential 

10  kmx  3  km 

1.0% 

3, 5, 7, 8, 9 

For  each  group  of  random  heterogeneity  parameters,  we  generate  5  velocity 
realizations  using  different  random  seeds,  and  calculate  the  wave  propagations  in  these 
models  for  up  to  1,000  km  distance.  The  source  time  function  is  a  Gaussian  derivative 


0-f)exp  -(*-*,)' 


,  with  t  =  0.2 sec.  The  earth  flattening  transform  (EFT)  is  used  to 


convert  the  spherical  earth  into  a  flat  earth.  A  factor  of  l/sqrt(r)  is  used  to  transfer  the 
geometry  from  2D  to  3D.  Illustrated  in  Figure  34  are  examples  of  Pn  waves  propagating  in 
models  without  and  with  random  velocity  fluctuations.  Figure  34a  and  34b  are  Pn  waves 
propagating  in  the  background  velocity  model,  and  34c  and  34d  are  Pn  wave  propagations 
in  a  velocity  model  with  1%  rms  velocity  perturbations.  The  pervasive  generation  of 
scattered  waves  is  readily  evident  in  the  snapshots,  particularly  in  the  Pn  coda. 


34 

Approved  for  public  release;  distribution  is  unlimited. 


distance  km 


distance  km 


distance  km 


distance  km 


Figure  34.  Wavefield  snapshots  in  velocity  models  without  (a),  (b)  and  with  (c),  (d) 
random  velocity  perturbations.  The  model  is  composed  of  a  40  km  thick  constant  velocity 
crust  and  an  upper  mantle  with  0.001  vertical  velocity  gradient.  The  upper  mantle  for  the 
calculations  in  (c)  and  (d)  has  1%  rms  velocity  perturbations  and  the  horizontal  and 
vertical  correlations  lengths  are  20  km  and  6  km,  respectively. 


The  synthetic  seismograms  are  processed  using  the  following  procedure 
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(1)  Bandpass  filter  the  seismograms  with  central  frequencies  0.4,  0.8,  1.5,  3.0,  6.0  and 
10.0  Hz. 

(2)  Group  velocity  windowing  between  8.2  km/s  and  7.6  km/s  is  used  for  sampling  the 
Pn  signal. 

(3)  The  rms  amplitude  is  calculated  within  the  group  velocity  window. 

(4)  All  amplitudes  are  normalized  by  the  spectra  from  the  corresponding  ID  velocity 
model  calculations  at  a  reference  distance  of  300  km. 

(5)  The  mean  value  from  five  realizations  of  the  random  structure  is  calculated  and 
used  as  the  result. 

Figures  35  to  43  present  calculations  of  the  Pn  wave  spectral  amplitudes  versus 
distance  for  random  velocity  models  number  1  to  9  listed  in  Table  4.  The  vertical 
coordinate  is  the  spectral  amplitude,  and  the  horizontal  coordinate  is  the  distance.  In  each 
figure,  panels  (a)  to  (f)  are  for  frequencies  0.4,  0.8,  1.5,  3.0,  6.0  and  10.0  Hz,  respectively. 
Light  colored  symbols  are  amplitudes  calculated  from  individual  random  realizations. 
Solid  squares  with  error  bars  are  mean  values  and  standard  deviations  from  all  five 
realizations.  The  solid  diamonds  are  from  the  background  ID  velocity  model.  As  a 
comparison,  shown  in  panel  (g)  are  decay  curves  for  different  frequencies  in  the  1-D 
background  model.  Shown  in  panel  (h)  are  related  decay  curves  in  models  with  random 
velocity  perturbations  in  the  upper  mantle.  All  amplitudes  are  normalized  by  the  amplitude 
in  the  background  model  at  a  reference  distance  of  300  km.  Thus  the  source  spectrum  is 
removed. 

For  each  model,  the  Pn  amplitudes  for  the  heterogeneous  structure  are  enhanced 
relative  to  the  expected  behavior  for  the  ID  reference  structure  underlying  the 
heterogeneity.  This  is  a  direct  manifestation  of  the  dominant  effect  of  scattering  of  strong 
energy  out  of  the  lid  waveguide  into  the  crustal  Pn  phase  as  depicted  in  Figure  33c  and 
captured  in  the  snapshot  in  Figure  34c, d.  The  effect  increases  with  frequency,  noting  that 
the  spectral  levels  are  relative  to  the  reference  distance  of  300  km,  where  the  structural 
effects  and  scattering  effects  are  relatively  minor.  The  random  models  have  variability  as 
expected,  but  the  average  behavior  over  the  five  realizations  is  quite  stable.  Some  of  the 
small-scale  structure  in  the  amplitude  curves  may  begin  to  average  out  with  a  greater 
number  of  realizations,  but  some  of  it  may  persist  as  a  function  of  the  specific  anisotropy 
of  the  heterogeneity.  While  smooth  curvature  of  the  patterns  and  the  increase  with 
distance  for  the  stronger  gradient  models  is  suppressed  in  the  scattering  calculations,  the 
overall  frequency  dependence  of  the  models  is  actually  enhanced  in  all  cases  relative  to 
the  ID  reference  case.  This  emphasizes  that  the  1Hz  calculations  of  Avants  et  al.  (201 1) 
capture  some  of  the  basic  effect  of  the  heterogeneous  structure,  but  the  frequency 
dependence  is  not  lost  by  the  scattering.  The  actual  small-scale  structure  on  any  given 
path  is  not  known,  but  the  assertion  by  some  researchers  that  it  will  homogenize  the 
signals  such  that  no  overall  frequency  dependence  of  the  Pn  spreading  would  be 
manifested  and  a  single  frequency-independent  power-law  decay  relationship  is  sufficient 
is  clearly  not  supported  by  these  calculations  (the  first  systematic  treatment  that  we  are 
aware  of  for  the  Pn  phase). 
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Figure  35.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.000_ran_20kmx6km_0.5%  (Model  No  1  in  Table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  frequencies  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID 
background  model.  Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes 
for  the  random  velocity  model. 
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Figure  36.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.000_ran_20kmx6km_1.0%  (Model  No  2  in  Table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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Figure  37.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.000_ran_20kmx6km_2.0%  (Model  No  3  in  Table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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Figure  38.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.001_ran_20kmx6km_0.5%  (Model  No  4  in  Table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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Figure  39.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.001_ran_20kmx6km_1.0%  (Model  No  5  in  table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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Figure  40.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.002_ran_20kmx6km_1.0%  (Model  No  7  in  table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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Figure  41.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.001_ran_40kmxl0km_0.5%  (Model  No  7  in  table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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Figure  42.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
mode  Grad._0.001_ran_40kmxl0km_1.0%  (model  No  8  in  table  1).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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Figure  43.  Pn  wave  spectral  amplitude  versus  distance  calculated  from  random  velocity 
model  Grad._0.001_ran_10kmx3km_1.0%  (Model  No  9  in  Table  4).  (a)  to  (f)  are  spectral 
amplitudes  at  0.4,  0.8,  1.5,  3.0,  6  and  10  Hz.  (g)  Amplitudes  in  the  ID  background  model. 
Different  symbols  are  for  different  frequencies,  (h)  Mean  amplitudes  for  the  random 
velocity  model. 
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For  a  velocity  model  with  zero  mantle  velocity  gradient  and  horizontal  and  vertical 
correlation  lengths  of  20  km  and  6  km,  Figure  44  compares  the  effects  of  changing  rms 
velocity  perturbations  on  the  Pn  wave  amplitudes.  In  this  figure,  panels  (a)  to  (f)  are  for 
frequencies  0.4,  0.8,  1.5,  3.0,  6.0  and  10  Hz.  The  open  triangles,  circles  and  squares  are  for 
random  models  with  rms  velocity  perturbations  of  0.5%,  1.0%  and  2.0%,  respectively.  As 
a  comparison,  the  solid  squares  are  for  the  background  velocity  model.  We  see  that  the 
scattering  process  causes  prominent  amplitude  increase.  The  amplitude  increase  is 
proportional  to  the  rms  velocity  perturbations  and  is  frequency  dependent,  with  high 
frequency  signals  being  more  sensitive  to  the  scattering.  Due  to  the  sphericity  of  the 
model,  the  amplitude  versus  distance  curves  in  the  background  model  curved  up  at  large 
distance.  The  scattering  effect  modifies  the  decay  curves  and  they  show  less  curvature. 
Even  though  the  result  is  calculated  from  five  realizations,  there  are  apparent  fluctuations 
in  these  decay  curves,  indicating  the  randomness  of  the  scattering  process. 

In  Figure  45,  comparisons  are  made  for  velocity  models  having  the  same  horizontal  and 
vertical  correlation  lengths  (20km  x  6km)  and  the  same  rms  velocity  perturbation  of  1%, 
but  different  upper  mantle  velocity  gradients  of  0.000,  0.001  and  0.002.  Panels  (a)  to  (f) 
are  for  spectral  amplitudes  at  0.4,  0.8,  1.5,  3.0,  6.0  and  10.0  Hz.  Triangles,  squares  and 
circles  are  for  models  with  mantle  velocity  gradients  0.000,  0.001  and  0.002,  respectively. 
The  open  symbols  are  for  random  velocity  models.  As  comparisons,  the  solid  symbols 
give  amplitudes  in  the  related  background  velocity  models.  A  prominent  phenomenon  is, 
that  in  the  presence  of  small-scale  heterogeneities  in  the  mantle,  the  scattering  effect  will 
dominate  the  Pn  spreading  function.  The  decay  curves  become  less  sensitive  to  the  upper 
mantle  velocity  gradient.  With  scattering,  the  decay  curves  tend  to  approach  the  “power- 
law”  behavior.  A  possible  explanation  is,  at  moderate  distances,  scattering  deflects  more 
energy  from  the  upper  mantle  into  the  crust,  raises  the  concave  part  of  the  decay  curve.  At 
larger  distances,  scattering  along  the  upper  mantle  channel  effectively  attenuates  the  Pn 
energy  offseting  the  upward  increase  of  the  decay  curve  at  large  distances.  The  combined 
effect  makes  a  more  uniform  rate  of  decay,  but  it  is  not  the  frequency-independent  case  as 
found  for  a  conical  head  wave. 

Figure  46  compares  Pn  wave  spectral  amplitudes  calculated  from  random  velocity  models 
having  the  same  upper  mantle  velocity  gradients  of  0.001,  the  same  rms  velocity 
perturbation  of  1%,  but  different  horizontal  and  vertical  correlation  lengths.  Panels  (a)  to 
(f)  are  for  spectral  amplitudes  at  0.4,  0.8,  1.5,  3.0,  6.0  and  10.0  Hz.  The  open  triangles, 
squares  and  circles  are  for  models  with  correlation  lengths  10  km  x  3  km,  20  km  x  6  km 
and  40  km  x  10  km,  respectively.  The  solid  squares  are  for  related  background  velocity 
model.  Within  the  investigated  frequency  band,  the  Pn  wave  is  weakly  dependent  on  the 
size  and  aspect  ratios  of  heterogeneities.  The  perturbations  with  horizontal  and  vertical 
correlation  lengths  of  10  km  and  3  km  cause  more  scattering. 
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Figure  44.  Comparison  between  Pn  wave  spectral  amplitudes  calculated  from  random 
velocity  models  No  1,2,  and  3  in  Table  4.  These  models  have  the  same  velocity  gradient 
and  random  parameters  except  the  rms  velocity  perturbations  are  0.5%,  1.0%  and  2.0%. 
(a)  to  (f)  are  for  spectral  amplitudes  at  individual  frequencies.  Open  triangles,  squares  and 
circles  are  for  models  with  0.5%,  1.0%  and  2.0%  velocity  perturbations.  The  solid  squares 
are  for  the  background  model. 
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Figure  45.  Comparison  between  Pn  wave  spectral  amplitudes  calculated  from  random 
velocity  model  No  2,  5,  and  6  in  Table  4.  These  models  have  the  same  horizontal  and 
vertical  correlation  lengths  20  km  and  6  km,  and  rms  velocity  perturbation  of  1%,  but  the 
velocity  gradients  are  0.000,  0.001  and  0.002.  Panels  (a)  to  (f)  are  for  spectral  amplitudes 
at  individual  frequencies.  Triangles,  squares  and  circles  are  for  models  with  velocity 
gradients  0.000,  0.001  and  0.002.  The  open  symbols  are  for  random  velocity  models.  The 
solid  circles  are  for  the  related  background  velocity  models. 
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Figure  46.  Comparison  between  Pn  wave  spectral  amplitudes  calculated  from  random 
velocity  model  No  5,  8,  and  9  in  Table  4.  These  models  have  the  same  upper  mantle 
velocity  gradients  of  0.001,  the  same  rms  velocity  perturbation  of  1%,  but  their  horizontal 
and  vertical  correlation  lengths  are  different,  (a)  to  (f)  are  for  spectral  amplitudes  at 
individual  frequencies.  The  open  triangles,  squares  and  circles  are  for  models  with 
correlation  lengths  10  km  x  3  km,  20  km  x  6  km  and  40  km  x  10  km.  The  solid  squares  are 
for  the  related  background  model. 
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4.6  Pn  Geometric  Spreading  Functions  for  Continental  and  Oceanic  Paths 

Pn  observations  from  the  North  Korean  test  site  and  for  earthquakes  in  northeast  China 
have  been  studied.  Appendix  1  reproduces  a  paper,  Zhao  et  al.  (2015),  published  in 
Journal  of  Geophysical  Research,  which  was  partially  supported  by  this  contract.  The  Pn 
paths  from  the  North  Korean  test  include  continental  and  oceanic  paths  that  are 
characterized  by  very  different  crust  and  upper  mantle  structures,  which  appear  to 
considerably  affect  the  Pn  wave  propagation.  In  this  section  we  present  additional 
modeling  of  2D  wave  propagation  along  the  corresponding  paths  to  supplement  the 
findings  given  in  Appendix  1.  There  are  numerous  factors  that  may  contribute  to  the 
observed  Pn  amplitudes,  including  the  crustal  thickness,  the  upper  mantle  velocity 
gradient,  P-wave  attenuation  in  the  uppermost  mantle,  and  scattering  along  the 
propagation  path,  as  explored  with  generic  models  in  the  previous  sections.  When  the 
source  and  station  are  located  in  different  environments,  the  transitional  zone  between  the 
continental  and  oceanic  structures  can  also  play  an  important  role  in  Pn-wave  propagation 
efficiency,  as  suggested  by  the  2D  simulations  in  section  4.4.  We  simulate  Pn  wave 
propagation  along  these  paths  for  nuclear  tests  located  in  North  Korea  and  compare  the 
model  calculations  with  the  observed  data. 

Shown  in  Figure  47  are  locations  of  the  North  Korea  nuclear  test  site  (NKTS)  and 
stations  providing  Pn  observations.  Pn  waves  observed  in  mainland  China  primarily 
traverse  continental  paths,  while  Pn  observations  from  stations  in  Japan  traverse  the  Japan 
Sea,  involving  primarily  oceanic  paths.  Given  that  the  explosion  sources  are  virtually 
isotropic  and  the  observations  span  a  large  azimuth  range,  Zhao  et  al.  (2015)  compared  the 
Pn  spectral  amplitudes  along  the  combined  continental  and  oceanic  paths  at  1000  km 
epicentral  distance  (recorded  by  those  stations  near  the  large  circle  in  Figure  47).  After 
removing  the  Pn-wave  source  excitation  functions,  the  spectral  amplitudes  are  presented 
in  Figure  48  as  a  function  of  azimuth  from  the  source,  with  different  symbols  indicating 
spectral  amplitudes  from  three  NKT  explosions.  The  black,  blue  and  red  colors  indicate 
0.8,  7.0  and  10.0  Hz  data.  Solid  circles  with  error  bars  indicate  their  mean  values  and 
standard  deviations  obtained  within  30°  azimuth  windows.  Prominent  differences  can  be 
observed  for  Pn  waves  along  different  paths.  For  azimuths  from  60°  to  180°,  the  Pn  signals 
crossing  the  oceanic  path  are  strongly  frequency  dependent,  with  extremely  low  high- 
frequency  amplitudes.  For  azimuths  between  230°  and  280°,  the  wave  paths  are  through 
continental  China,  and  the  low-frequency  content  (0.8  and  7  Hz)  is  similar  to  that  along 
the  oceanic  paths,  but  the  high-frequency  (10  Hz)  propagation  is  much  more  efficient  than 
for  the  oceanic  paths.  The  geometric  spreading  function  tends  to  raise  the  high-frequency 
signal,  whereas  attenuation  tends  to  reduce  the  high-frequency  content.  The  detailed 
structure  of  transition  zones  may  also  affect  the  frequency  dependence. 


50 

Approved  for  public  release;  distribution  is  unlimited. 


50 


Figure  47.  Map  depicting  locations  of  the  North  Korean  test  site  (red  star  labeled  NKTS) 
and  seismic  stations  at  which  Pn  phase  measurements  were  made.  The  stations  located 
near  the  large  circle  have  nominal  epicentral  distance  of  1000  km  from  the  NKTS,  with 
the  pink  and  blue  colors  indicating  the  Pn  waves  that  traverse  primarily  continental  or 
oceanic  paths,  respectively.  (Zhao  et  al.,  2015  -  Appendix  1). 

To  simulate  the  propagation  of  Pn  waves  in  oceanic  and  continental  structures,  we 
designed  a  Japan  Sea  model  as  shown  in  the  top  panel  in  Figure  49.  The  depth  of  the 
Moho  discontinuity  is  obtained  from  CRUST1.0  (Laske,  et  al.,  2013)  by  using  a  typical  Pn 
path.  In  this  model,  the  crust  beneath  the  NKTS  is  approximately  33  km,  and  under  the 
Japan  Sea  is  approximately  11  to  15  km.  On  the  NKTS  side,  the  crustal  thickness  quickly 
decreases  from  33  km  to  about  10-12  km  over  approximately  100  km  distance,  forming  a 
steep  reduction  in  depth  of  the  Moho  discontinuity.  On  the  side  of  the  Japan  Islands,  the 
crustal  thickness  increases  relatively  slowly,  forming  a  broad  dip  in  the  Moho 
discontinuity.  For  continental  path,  we  use  a  simple  1-D  model  with  a  30  km  flat  Moho 
(bottom  panel  in  Figure  49).  In  both  models,  the  upper  mantle  velocity  gradient  is  set  to  be 
zero,  and  the  water  column  is  not  included  in  the  model. 
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Figure  48.  Pn  spectral  amplitudes  measured  from  stations  near  a  distance  of  1000  km 
versus  azimuth  (Figure  47).  Different  symbols  are  amplitudes  from  the  three  North  Korea 
nuclear  tests.  Black,  blue  and  red  colors  indicate  0.8,  7.0  and  10.0  Hz  measurements.  Solid 
circles  with  error  bars  represent  their  mean  values  and  standard  deviations,  obtained 
within  30-degree  azimuth  windows.  The  Pn  source  excitation  functions  are  removed  from 
the  data.  (Zhao  et  al.,  2015  -  Appendix  1). 
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Figure  49.  Crust-upper  mantle  models  for  Pn-wave  simulation.  Top:  a  typical  profile  from 
the  NKTS  to  Japan  Island  passing  the  Japan  Sea.  Bottom:  a  30  km  thick  crust  model  for 
typical  continental  path.  To  show  the  details,  the  vertical  coordinate  is  exaggerated. 
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To  calculate  the  synthetic  seismograms,  we  use  the  same  source  time  function  as  used  in 
the  previous  calculations.  However,  we  set  the  source  depth  at  1  km.  Shown  in  Figure  50 
are  wavefield  snapshots  for  the  entrance  section  near  the  NKTS  (from  50  km  to  350  km) 
of  the  Japan  Sea  model.  Shown  in  Figure  5 1  are  wavefield  snapshots  for  the  Japan  Island 
side  (from  650  km  to  950  km).  Shown  in  Figures  52  and  53  are  similar  snapshots  for  Pn 
waves  in  the  continental  model  which  has  a  30  km  thick  crust.  A  shallowing  Moho  tends 
to  allow  more  energy  from  the  source  to  penetrate  into  the  upper  mantle,  while  a 
deepening  Moho  allows  more  energy  to  come  back  from  the  upper  mantle  into  the  crust. 
Aside  from  any  intrinsic  attenuation,  these  conditions  favor  Pn  wave  traversing  the  Japan 
Sea  model  compared  to  transmission  along  a  flat  continental  path.  This  can  also  be 
understood  by  comparing  velocity  models  in  Figure  49,  where  the  varying  dip  of  the  Japan 
Sea  Moho  is  equivalent  to  an  additional  curvature  relative  to  the  1-D  model.  Figure  54 
compares  the  Pn  waves  in  both  models  when  they  approach  1000  km  distance.  It  is 
evident  that  in  the  Japan  Sea  model  Pn  is  stronger  than  in  the  continental  model. 
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Figure  50.  Wavefield  snapshots  of  Pn  waves  propagating  in  the  entrance  section  near  the 
NKTS  (50  -  350  km)  of  the  Japan  Sea  model. 
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Figure  51.  Wavefield  snapshots  of  Pn  waves  propagating  in  the  exit  section  approaching 
the  Japan  receivers  (650  -  950  km)  of  the  Japan  Sea  model. 
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Figure  52.  Wavefield  snapshots  of  Pn  waves  propagating  in  the  entrance  section  (50  -  350 
km)  of  the  1-D  model  with  a  30  km  thick  crust. 
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Figure  53.  Wavefield  snapshots  of  Pn  waves  propagating  in  the  exit  section  (650  -  950 
km)  of  the  30  km  thick  crust  model. 
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Figure  54.  Comparison  between  Pn  waves  propagating  in  the  continental  and  Japan  Sea 
models  between  650  and  950  km.  Note  that  the  Pn  wave  amplitude  is  stronger  in  the  Japan 
Sea  model. 
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Figure  55  shows  the  synthetic  seismograms  for  both  the  Japan  Sea  model  and  the  flat 
continental  model,  where  from  left  to  right  are  the  Japan  Sea  velocity  model,  synthetics 
for  the  Japan  Sea  model  and  synthetics  for  the  continental  model.  The  Pn  wave  amplitude 
is  much  stronger  in  the  Japan  Sea  model.  Because  the  Pn  wave  penetrates  the  low-velocity 
crust  with  laterally  varying  thickness,  its  travel  time  curve  deviates  from  a  straight  line. 
The  vertical  markers  on  the  seismograms  indicate  the  group  velocity  windows  used  to 
measure  the  rms  amplitude.  The  frequency  dependent  Pn  amplitudes  versus  distance  are 
illustrated  in  Figure  56,  where  the  diamonds  are  for  the  Japan  Sea  model,  crosses  are  for  a 
30  km  thick  flat  crust  model,  and  plus  signs  are  for  a  10  km  thick  flat  crust  model.  It  is 
evident  that  the  crustal  thickness  does  not  affect  the  Pn  amplitude  much.  However,  the  Pn 
amplitudes  in  the  Japan  Sea  model  are  stronger  than  those  in  flat  Moho  models.  These 
differences  are  frequency  dependent.  At  lower  frequencies,  the  difference  is  less  than  a 
half  order  of  magnitude,  but  it  increases  to  about  an  order  of  magnitude  at  higher 
frequencies.  It  appears  the  deepening  Moho  on  the  Japan  Island  side  also  plays  an 
important  role  in  changing  the  Pn  wave  amplitudes.  At  lower  frequencies  it  raises  the  Pn 
amplitudes.  For  higher  frequencies,  its  effect  is  complicated. 


Japan  Sea  model 


30  km  crust 


Figure  55.  Comparison  between  synthetic  seismograms  in  the  Japan  Sea  model  and  a  1-D 
continental  model  with  a  30  km  thick  crust.  Left:  Japan  Sea  velocity  model,  middle: 
synthetic  seismograms  in  Japan  Sea  model,  and  right:  synthetic  seismograms  in  the  model 
with  a  30  km  thick  crust.  The  vertical  markers  indicate  the  group  velocity  windows  used 
to  measure  the  rms  amplitude. 
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Figure  56.  Frequency  dependent  amplitude  versus  distance  for  the  Japan  Sea  model 
(diamond),  30  km  crust  (cross)  and  10  km  crust  (plus  sign)  models. 


Given  that  the  oceanic  dataset  is  mostly  observed  in  Japan  over  a  narrow  distance  range  at 
approximately  1000  km  distance,  we  cannot  compare  Pn  amplitude  variations  versus 
distance  with  our  simulations.  Instead,  we  investigate  the  observed  frequency  dependence 
of  Pn  wave  amplitudes  at  1000  km.  Similar  to  the  measurements  shown  in  Figure  48,  we 
A  ( f) 

calculate  the  Pn-mokm\J  )  from  the  NKTS  dataset.  Figure  57  shows  spectral  amplitudes 
for  both  the  oceanic  path  through  the  Japan  Sea  and  the  continental  path  through  mainland 
China.  At  lower  frequencies,  Pn  waves  from  both  paths  have  similar  amplitudes. 
However,  at  higher  frequencies,  the  Pn  waves  traversing  the  Japan  Sea  path  are  much 
weaker  than  on  continental  paths  towards  mainland  China.  We  know  that  the  geometric 
spreading  tends  to  raise  the  high  frequency  signal  while  the  attenuation  tends  to  decrease 
the  high  frequency  signal.  By  comparing  the  numerical  simulations  of  Pn  wave 
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geometrical  spreading  (Figures  55  and  56)  and  the  observations  from  the  NKTS  data 
(Figure  57),  we  see  that,  even  though  the  geometric  spreading  in  the  Japan  Sea  model 
appears  very  efficient,  the  observed  Pn  amplitudes  at  high  frequency  through  the  Japan 
Sea  is  weaker.  This  suggests  that  the  P  wave  propagation  in  the  uppermost  mantle  beneath 
the  Japan  Sea  must  be  more  attenuating  than  that  under  the  continental  crust.  To  compare 
with  the  observed  Pn  wave  spectra,  we  combine  Pn  Q  models  with  the  calculated 
geometrical  spreading  functions.  For  Pn  waves  from  the  NKTS  to  continental  stations, 

O  =237  f036 

Zhao  et  al.  (2015)  obtained  a  model  of  ^ Pn  ’  ,  which  will  be  used  for  the 

continental  path.  For  the  oceanic  model,  using  a  trial  and  error  method,  we  find  a  lower 

O  O  =150/’° 40 

and  a  slightly  higher  r|  giving  a  Pn  Q  model  of  ^Pn  J  for  the  Japan  Sea.  The 

resulting  Pn  wave  amplitudes  are  shown  in  Figure  58.  The  predictions  are  consistent  with 

the  observed  data  shown  in  Figure  57.  As  discussed  above,  there  are  many  factors 

affecting  the  Pn  wave  amplitude  and  its  frequency  dependence.  Thus,  this  result  is  not 

unique.  However,  it  provides  some  constraint  on  the  P  wave  attenuation  beneath  the 

oceanic  crust. 


Figure  57.  Observed  Pn  wave  amplitudes,  at  1000  km  distance  from  the  NKTS,  versus 
frequency.  Solid  diamonds  are  for  paths  traversing  the  Japan  Sea,  and  open  circles  are  for 
paths  traversing  continental  structure  to  Northwest  China.  Pn  waves  through  these  two 
paths  show  apparently  different  frequency  dependence. 
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Figure  58.  Computed  Pn  wave  amplitudes  at  1000  km  distance  from  the  NKTS  versus 
frequency.  Solid  diamonds  are  for  the  Japan  Sea  model,  and  open  circles  are  calculated  in 
the  30  km  thick  crustal  model  for  the  continental  path.  The  results  are  obtained  by 
combining  the  calculated  geometrical  spreading  with  Pn  Q  models.  For  the  continental 
model,  the  Pn  Q  model  QPn  =237 /°'36  (Zhao  et  al.,  2015)  is  used.  For  the  Japan  Sea 

model,  a  Pn  Q  model  of  QPn  =  150 /040  is  used. 

4.7  Use  of  RSTT  Models  and  Empirical  Models  for  Pn  Spreading  for  Eurasian  Data 

In  working  toward  self-consistent  use  of  velocity  models  for  predicting  regional  phase 
travel  times  and  regional  phase  geometric  spreading,  we  empirically  examined  the  use  of 
regionalized  corrections  for  Pn  amplitudes  guided  by  RSTT  regional  velocity  models  and 
by  empirical  analysis  of  Pn  data  with  constraint  of  realistic  spreading  behavior  as 
established  by  Yang  (2011).  Figure  59  displays  the  regionalization  guided  by  RSTT 
models.  It  would  be  a  major  computational  effort  to  compute  the  fully  3D  models  for  each 
region  to  evaluate  Pn  spreading,  and  the  results  of  earlier  sections  have  established  that 
laterally  varying  gradients  and  crustal  thickness  produce  directionally-dependent  path 
corrections,  so  comparisons  with  data  are  very  complex  and  essentially  one  would  have  do 
to  a  path-by-path  treatment.  This  is  unwieldy,  and  we  adopt  a  simplified  procedure  for 
evaluating  whether  such  an  approach  is  even  warranted.  The  RSTT  models  are  obtained 
from  noisy  and  limited  data  sets  that  provide  only  limited  constraint  on  mantle  velocity 
gradients  (sometimes  they  may  capture  true  travel  time  curvature  and  constrain  the 
gradient,  but  in  most  cases  the  data  are  too  noisy  and/or  sparse  to  really  resolve  actual  lid 
velocity  gradient.  Our  approach  is  to  explore  first-order  variations  in  structure  inferred 
from  the  RSTT  model  variations. 
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Figure  59.  RSTT  model  lid  velocity  gradients  (top)  and  average  Moho  P  velocity  (bottom) 
with  a  superimposed  grid  developed  to  define  provinces  of  similar  parameters. 

The  LANL  database  for  Pn  measurements  (Figure  60)  is  used  to  extract  signals  and 
measured  Pn  amplitudes  in  varying  passbands  to  define  the  sampling  of  regions  (Figure 
61)  for  which  RSTT  mantle  lid  gradient  and/or  average  Moho  velocity  are  relatively 
uniform,  or  at  least  we  approximate  with  a  regional  ID  structure. 
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Figure  60.  LANL  Eurasian  Pn  data  base  coverage  showing  Pn  paths  sampling  the  RSTT- 
based  regionalization. 
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Figure  61.  Regional  numbering  scheme  for  the  RSTT -based  subdivision.  Pn  observations 
in  each  region  have  been  evaluated.  Model  spreading  parameters  for  each  RSTT  structure 
have  been  determined. 


For  each  region  identified  in  Figure  61,  we  computed  spreading  behavior  for  a  ID 
average  of  the  local  RSTT  model  using  wavenumber  integration,  and  applied 
corresponding  frequency-dependent  geometric  spreading  corrections  to  the  MDAC- 
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corrected  Pn  observations  for  that  region.  We  also  applied  simple  power-law  spreading 
corrections  for  a  power  of  -1.1,  and  we  determined  region-specific  spreading  corrections 
following  the  procedure  of  Yang  (2011).  This  gave  three  sets  of  corrected  Pn  data  for  each 
region.  Plotting  these  sets  of  corrected  Pn  amplitudes  (log  10)  as  functions  of  distance, 
allows  an  estimation  of  the  average  Q  for  that  region  for  each  frequency.  We  display 
results  for  regions  4  and  13  (Figure  61)  in  Figures  63a, b  and  64a, b,  respectively. 
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Figure  62a.  Pn  amplitudes  versus  distance  for  Region  4  with  corrections  using  a  power- 
law  (top  row),  RSTT-based  synthetic  computations  (middle  row)  or  empirical  region- 
specific  corrections  following  the  procedure  of  Yang  (2011)  (bottom  row)  for  passband 
center  frequencies  of  1.06  Hz  (left)  and  4.24  Hz  (right).  The  measured  average  Q  for  each 
case  obtained  from  the  regression  curve  is  indicated  in  the  upper  right. 

Each  set  of  geometric  spreading  corrections  results  in  relatively  smooth  linear 
amplitude  versus  distance  trend  that  suggests  a  corresponding  average  Q  for  that 
frequency.  However,  the  Q  values  that  are  inferred  vary  dramatically,  depending  on  the 
geometric  spreading  used. 
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Region  4:  Power-law  spreading  correction 


Region  4:  Region-specific  spreading  correction 
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Figure  62b.  Pn  amplitudes  versus  distance  for  Region  4  (map  inset)  with  corrections  using 
a  power-law  (top  row),  RSTT-based  synthetic  computations  (middle  row)  or  empirical 
region-specific  corrections  following  the  procedure  of  Yang  (2011)  (bottom  row)  for 
passband  center  frequency  of  8.94  Hz  (right).  The  measured  average  Q  for  each  case 
obtained  from  the  regression  curve  is  indicated  in  the  upper  right. 

In  general,  as  found  by  Yang  et  al.  (2007)  and  Yang  (2011),  use  of  a  ‘standard’  power-law 
spreading  correction  with  an  exponent  of  -1.1  that  does  not  vary  with  frequency  (and  for 
which  there  is  no  specific  known  structure)  gives  high  Q  values  in  every  case,  and  even 
negative  Q  value  for  Region  13  at  1.06  Hz.  For  Region  4,  the  RSTT  velocity  gradient  is 
relatively  high  (0.0026  s’1)  while  for  Region  13  it  is  average  value  (0.0012  s'1).  In  both 
cases,  the  addition  of  the  EFT  gradient  of  0.0013  s'1  gives  positive  velocity  gradients  that 
produce  strong  frequency-dependent  Pn  amplitude  corrections.  The  data  have  somewhat 
less  scatter  and  systematically  lower  Q  values  are  inferred  for  each  region  using  the  RSTT 
based  corrections.  Q  values  that  are  1.6- 1.8  times  higher  than  the  RSTT-based  results  are 
found  for  corrections  based  on  the  region-specific  spreading  corrections  following  Yang 
(201 1).  The  absolute  values  and  the  frequency-dependence  of  the  attenuation  parameters 
are  more  similar  for  the  RSTT-  and  regional-based  spreading  corrections.  For  these 
regions  there  is  no  clear  reduction  in  scatter  to  favor  a  specific  spreading  model,  but  the  Q 


64 

Approved  for  public  release;  distribution  is  unlimited. 


values  appear  to  be  more  reasonable  for  either  of  the  frequency-dependent  spreading 
corrections  results.  As  long  as  the  Q  model  for  a  given  spreading  model  is  used,  total 
amplitude  calculations  will  be  reasonably  similar.  The  main  advantages  of  using  a  more 
appropriate  local  velocity  model  to  process  the  amplitudes  is  that  a  more  reliable  Q  model 
may  be  estimated  and  that  could  be  valuable  for  application  to  full  waveform  modeling, 
for  extrapolation  to  poorly  sampled  paths,  and  for  overall  self-consistency  of  event 
locations  and  amplitudes. 
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Figure  63a.  Pn  amplitudes  versus  distance  for  Region  13  with  corrections  using  a  power- 
law  (top  row),  RSTT-based  synthetic  computations  (middle  row)  or  empirical  region- 
specific  corrections  following  the  procedure  of  Yang  (2011)  (bottom  row)  for  passband 
center  frequencies  of  1.06  Hz  (left)  and  4.24  Hz  (right).  The  measured  average  Q  for  each 
case  obtained  from  the  regression  curve  is  indicated  in  the  upper  right. 
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Region  13:  Power-law  spreading  correction 


Figure  63b.  Pn  amplitudes  versus  distance  for  Region  13  (map  inset)  with  corrections 
using  a  power-law  (top  row),  RSTT-based  synthetic  computations  (middle  row)  or 
empirical  region-specific  corrections  following  the  procedure  of  Yang  (2011)  (bottom 
row)  for  passband  center  frequency  of  8.94  Hz  (right).  The  measured  average  Q  for  each 
case  obtained  from  the  regression  curve  is  indicated  in  the  upper  right. 

These  comparisons  support  the  premise  that  using  regionally  varying  velocity  models 
can  enable  self-consistent  travel  time  and  geometric  spreading  for  the  same  regional 
phases.  The  scatter  is  large  though,  and  one  cannot  demonstrate  that  a  specific  regional 
model  and  its  associated  prediction  of  geometric  spreading  performs  most  reliably  in 
terms  of  reducing  variance.  This  is  not  unexpected  given  the  results  of  the  earlier  sections 
demonstrating  that  in  the  presence  of  laterally  varying  lid  structure  the  actual  path 
direction  in  the  medium  has  a  significant  effect  on  both  the  absolute  spreading  correction 
and  the  frequency-dependent  trends  in  the  spreading  with  distance.  It  is  pretty  clear  that  if 
the  structure  is  well-determined  by  travel  time  branch  curvature  and/or  modeling  of 
waveforms,  use  of  that  structure  for  regional  phase  spreading  is  the  most  sensible 
approach  (compared  to  a  frequency-independent  power  law  correction  with  an  attendant, 
biased  frequency-dependent  attenuation  model).  But,  the  reality  is  that  neither  RSTT  nor 
other  regional  studies  provide  strong  constraints  on  actual  lid  velocity  models  at  this  time. 


66 

Approved  for  public  release;  distribution  is  unlimited. 


4.8  A  2D  Pn  Velocity  Gradient  Model  for  Eurasia 

P-wave  velocities  in  the  uppermost  mantle  change  with  depth.  The  rate  of  this  velocity 
change,  or  the  radial  P- wave  velocity  gradient,  plays  an  important  role  in  shaping  the  Pn 
propagation  behavior  both  in  terms  of  its  travel  time  and  its  amplitude  as  demonstrated  in 
the  previous  sections.  The  velocity  gradient  also  varies  laterally  due  to  different 
uppermost-mantle  conditions.  It  may  reflect  lateral  variations  of  temperature,  composition, 
and  dynamics  of  the  uppermost  mantle. 

Zhao  (1993)  mathematically  formulated  P„  travel  time  as  the  summation  of  crustal-leg 
travel  time,  head-wave  travel  time  and  the  contribution  from  Pn  velocity  gradient.  Using 
this  formulation,  Zhao  (1993)  and  Zhao  and  Xie  (1993)  estimated  the  average  Pn  velocity 
gradient  for  the  Basin  and  Range  province  in  the  United  States  and  for  the  Tibetan  Plateau 
in  China.  Phillips  et  al.  (2007)  and  Myers  et  al.  (2010)  extended  Zhao’s  (1993) 
formulation  to  account  for  the  lateral  variation  of  P„  velocity  gradient  and  created  two- 
dimensional  (2D)  P„  velocity-gradient  models  for  Eurasia  and  North  Africa  that  are 
merged  with  locally  determined  regional  velocity  structures  in  the  RSTT  formulation. 

In  this  section,  we  adopted  the  formulation  of  Zhao  (1993)  to  estimate  path-specific  Pn 
velocity  gradients  from  observed  P„  travel  times.  We  assumed  that  each  path-specific 
velocity  gradient  is  the  mean  of  laterally  varying  velocity  gradients  along  the  path.  Using 
this  assumption,  we  inverted  the  path-specific  P„  velocity  gradients  for  a  2D  Pn  velocity- 
gradient  model  for  Eurasia.  This  model  framework  is  helpful  for  extracting  models  to  be 
used  in  Pn  geometric  spreading  calculations  for  specific  paths,  as  an  alternative  to  the 
subregion  averaging  of  RSTT  models  used  in  the  previous  section.  Our  development  of  a 
2D  Pn  velocity  gradient  model  for  Eurasia  utilizes  different  path  averaging  that  was  used 
in  RSTT,  and  we  show  that  our  approach  may  provide  a  better  averaging  of  the  data 
behavior,  although  the  models  are  still  data-limited  with  respect  to  robust  determination  of 
radial  velocity  structure. 

We  selected  data  from  an  updated  P„  travel-time  dataset  of  Myers  et  al.  (2010) 
(Begnaud,  person,  comm.).  It  contains  Pn  travel  times  from  more  than  5,000  earthquakes 
and  explosions  recorded  by  about  4,500  stations  throughout  Eurasia.  Epicentral  distances 
of  the  dataset  are  between  2°  and  15°.  We  required  that  events  for  which  we  include  Pn 
data  should  have  estimated  depths  (not  fixed  depths)  and  they  should  not  be  in  the  mantle 
(depth  <  40  km).  The  total  number  of  selected  measurements  is  about  85,000. 

Figure  64  plots  the  path  density  of  the  dataset  that  we  used  to  invert  for  the  velocity- 
gradient  model.  In  general,  P„  path  density  is  high  in  tectonically  active  regions  due  to 
large  numbers  of  earthquakes  and  stations  in  those  regions.  This  is  particularly  true  for 
southern  Europe,  where  the  P„  path  density  is  the  highest.  The  large  number  of  paths 
provides  relatively  uniform,  albeit  sparse  sampling  of  all  of  Eurasia  as  indicated  in  Figure 
64. 
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Figure  64.  Path  density  of  Pn  travel-time  measurements  used  in  the  tomography. 

Zhao  (1993)  formulated  P„  travel  time  from  event  i  to  station  j  as 

tg=ti+tj  +  th+y  ’ 

where  ty  is  the  total  Pn  travel  time,  t,  is  the  vertical  Pn  travel  time  in  the  crust  below  the 
source  and  tj  is  the  vertical  P„  travel  time  in  the  crust  below  the  station.  Zhao  (1993) 
decomposed  the  rest  of  the  Pn  travel  time  into  4  and  y,  where  4  was  termed  “horizontal 
travel  time”  by  Zhao  (1993)  and  is  composed  of  horizontal  Pn  travel  times  in  the  crust 
below  the  source  and  the  station,  and  the  Pn  travel  time  along  the  Moho  between  Moho 
piercing  points  of  the  Pn  path,  y  is  the  contribution  of  Pn  velocity  gradient  to  the  travel 
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time.  Myers  et  al.  (2010)  modified  Eq.  (4)  and  redistributed  P„  horizontal  travel  times  in 
the  crust  from  4  to  t,  and  tj  such  that  4  contained  only  P„  horizontal  travel  time  along  the 
Moho  between  Pn  piercing  points.  Under  the  assumption  that  P„  velocity  gradient  was 
constant  for  a  particular  region  and  the  gradient  g  was  much  smaller  than  the  average  P- 
wave  velocity  vo  beneath  the  Moho,  Zhao  (1993)  derived  the  following  expression  for  y 


7 


c2X3 
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where  c  =  g/v 0  and  X  is  the  horizontal  distance  between  Pn  Moho  piercing  points.  For  a 
spherical  Earth,  c  =  (g/v 0  +  1/r)  where  r  is  the  distance  from  the  center  of  the  Earth  to 
Moho.  1/r  is  an  Earth  flattening  correction. 


To  obtain  a  2D  Pn  velocity- gradient  map  as  well  as  a  2D  Pn  velocity  map  using  Eq. 
(4),  Phillips  et  al.  (2007)  and  Myers  et  al.  (2010)  assumed  that  c2  in  Eq.  (5)  could  be 
represented  as 
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where  c2  are  discretized,  laterally- varying  c2  values  for  the  region  and  x„  is  the  length  of 

the  nth  path  segment  .  As  a  result,  Eq.  (4)  was  linearized  and  a  linear-inversion  method 
could  be  used  for  the  tomography. 

In  this  study,  we  made  a  different  assumption.  Instead  of  representing  the  path-specific 
c2  as  the  mean  of  laterally- varying  c2,  we  assumed  that  the  path-specific  velocity  gradient 
g  itself  was  the  mean  of  the  laterally-varying  velocity  gradient  along  the  Pn  path.  That  is, 

S  =  (y) 

n=l 

To  make  use  of  Eqs.  (4)  and  (7)  in  a  linear  inversion  to  obtain  a  2D  Pn  velocity- 
gradient  model,  we  need  to  subtract  the  first  three  terms  on  the  right-hand  side  of  Eq.  (4) 
from  the  left-hand  side.  The  velocity  gradient  g  for  a  particular  path  can  then  be  calculated 
with  Eq.  (5)  and  used  in  Eq.  (7)  for  the  inversion.  We  achieved  this  by  tracing  the  Pn  path 
through  the  Regional  Seismic  Travel  Time  (RSTT)  model  (Myers  et  al.,  2010)  and 
calculating  the  travel  times  4  tj  and  4.  Because  RSTT  was  obtained  from  an  inversion  that 
incorporated  both  Pn  slowness  and  Pn  velocity  gradient  as  model  parameters,  and  these 
two  parameters  depend  on  the  path  in  similar  ways  (Myers  et  al.,  2010,  Eq.  7),  there  might 
be  certain  degree  of  tradeoffs  between  the  resulting  velocity  map  and  the  velocity-gradient 
map.  To  avoid  the  potential  error  of  the  RSTT  Pn  velocity  map  due  to  the  tradeoff,  we 
replaced  the  Moho  Pn  velocity  in  RSTT  with  that  from  the  independent  Sandia-Los 
Alamos  3D  P- wave  velocity  model  (SALSA3D)  (Ballard  et  al.,  2010;  Begnaud  et  al., 
2011)  and  constructed  a  hybrid  model  for  the  ray  tracing. 
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Figure  65  shows  the  calculated  y  by  subtracting  U,  tj  and  4,  obtained  from  the  ray 
tracing,  from  measured  Pn  travel  times  (Eq.  4).  According  to  Eq.  (5),  y  should  always  be 
negative  and  decreases  with  increasing  distance.  A  smaller  y  corresponds  to  a  higher 
velocity  gradient  in  the  uppermost  mantle. 
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Figure  65.  Path-specific  y  from  measured  P„  travel  times. 

To  see  which  assumption,  either  the  Phillips  et  al.  (2007)  and  Myers  et  al.  (2010) 
assumption  of  path-specific  c2  as  the  mean  of  location-dependent  c2  or  our  assumption  of 
path-specific  gradient  as  the  mean  of  location-dependent  gradient,  is  more  reasonable,  we 
conducted  two-dimensional  Monte  Carlo  simulations  to  simulate  Pn  waves  traveling  in  a 
mantle  with  laterally  heterogeneous  velocity  gradient. 

The  2D  model  that  we  used  for  the  simulation  was  composed  of  10  equal- width 
regions  with  the  same  velocity  of  8.1  km/s  at  the  top  and  different  vertical  velocity 
gradients.  The  model  is  500  km  in  depth,  which  guarantees  that  all  Pn  rays  bottom  within 
the  model.  In  the  simulation,  travel  times  between  a  source  and  a  receiver  1200-km  apart 
at  the  two  top  comers  of  the  model  were  calculated  repeatedly.  For  each  calculation,  the 
order  of  the  regions  was  changed  in  a  random  fashion  from  the  previous  calculation. 
Figure  66  shows  two  realizations  of  the  model  as  examples. 
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Figure  66.  Two  realizations  of  the  velocity-gradient  model  for  the  Monte  Carlo 
simulation. 

We  used  the  Multistencils  Fast  Marching  (MSFM)  method  (Hassouna  and  Farag, 
2007)  to  calculate  Pn  travel  times  through  the  model.  Figure  67  plots  the  P„  travel  times 
calculated  from  10,000  iterations  as  a  histogram.  The  red  line  in  the  figure  marks  the 
travel  time  calculated  from  a  model  with  a  single  velocity  gradient  that  is  the  mean  of  the 
gradients  of  the  model  used  in  the  simulation.  The  dashed  green  line  is  the  travel  time 
from  a  model  constructed  using  the  assumption  of  Phillips  et  al.  (2007)  and  Myers  et  al. 
(2010).  The  travel  time  calculated  using  our  assumption  is  very  close  to  the  mean  of  the 
travel  times  from  the  Monte  Carlo  simulation  (0.07%  difference).  The  travel  time 
calculated  from  Phillips  et  al.  (2007)  and  Myers  et  al .’  (2010)  assumption,  on  the  other 
hand,  is  much  more  biased  (0.7%  difference).  The  simulation  indicates  that  the  cumulative 
velocity-gradient  effect  on  the  Pn  travel  time  is  better  modeled  by  averaging  the  laterally- 
varying  velocity  gradient  along  the  Pn  path  than  by  averaging  the  square  of  the  gradient. 
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Figure  67.  Distribution  of  Pn  travel  times  from  the  Monte  Carlo  simulation.  The  mean  of 
the  distribution  is  143.4s.  The  red  line  marks  the  travel  time  for  a  model  with  a  single 
velocity  gradient  calculated  as  the  mean  of  the  laterally  varying  gradients.  The  green  line 
marks  the  travel  time  for  a  model  where  the  gradient  is  calculated  as  the  square  root  of  the 
mean  of  the  square  of  the  laterally  varying  gradients. 

We  took  a  Bayesian  approach  in  a  damped  and  weighted  least-squares  inversion  to 
invert  for  the  2D  Pn  velocity- gradient  model.  We  discretized  Eurasia  using  the  software 
GeoTess  (Ballard  et  al.,  2009;  Ballard  et  al.,  2012),  which  employs  a  triangular 
tessellation  method.  For  such  a  node-based  model,  Eq.  (7)  becomes 

Y  N  3 
A  u=l  p=i 

where  gp  are  velocity  gradients  at  model  nodes  surrounding  the  nth  path  segment  and  cpn 
are  weighting  coefficients  based  on  the  distances  from  the  center  of  the  path  segment  to 
the  surrounding  model  nodes.  With  calculated  path-specific  P„  velocity  gradients  and  Eq. 
(8),  we  built  a  system  of  linear  equations  for  the  tomographic  inversion  to  construct  the 
2D  velocity-gradient  model. 

Figure  68  illustrates  the  model  tessellation.  To  construct  the  tessellation,  we  first 
defined  the  region  where  we  had  data  coverage.  We  then  constructed  the  mesh  using 
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variable-resolution  discretization.  Within  the  region  with  data  coverage,  the  node  spacing 
is  1°.  Away  from  the  region,  the  spacing  progressively  increases.  During  the  inversion, 
only  nodes  that  contributed  to  the  calculation  of  path-specific  velocity  gradients  were 
included. 


Figure  68.  Tessellation  of  the  model. 

We  used  the  right-hand  side  of  Eq.  (5)  to  derive  an  a  priori  model  for  the  inversion. 
We  first  summed  cpn  for  each  node  over  all  paths.  We  then  replaced  gp  on  the  right-hand 
side  of  the  equation  with  g  on  the  left-hand  side  for  each  measurement  and  summed  gcpn 
for  each  node  over  all  paths.  The  a  priori  model  was  obtained  by  dividing  the  second 
summation  by  the  first  summation.  The  a  priori  model  is  shown  in  the  Figure  69.  The 
gradients  range  from  -0.001  s'1  to  0.074  s'1. 
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Figure  69.  The  a  priori  model  used  in  the  inversion. 

Figure  70  shows  the  P„  velocity-gradient  model  for  Eurasia  from  the  tomographic 
inversion  along  with  constant  model-resolution  contours.  A  model  node  is  perfectly 
resolved  if  its  resolution  is  one.  The  best-resolved  node  in  the  model  has  a  resolution  of 
0.73.  The  model  is  reasonably  resolved  along  the  Alpine-Himalayan  orogenic  belt,  in 
Scandinavia,  and  along  the  western  Pacific  in  Japan  and  Taiwan.  Northern  Eurasia, 
Arabian  Peninsula,  southeast  China  and  much  of  India  are  not  very  well  resolved  due  to 
poor  path  coverage. 
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Figure  70.  P„  velocity-gradient  model  from  the  tomographic  inversion.  Purple  lines  are 
constant  model-resolution  contours.  Corresponding  values  of  some  of  the  contours  are 
marked. 

To  better  visualize  the  resolution  power  of  the  data  used  in  the  inversion,  Figure  71 
shows  the  result  of  a  checkerboard  test.  A  checkerboard  model  with  5°  cells  and 
alternating  gradient  values  between  -0.0006  s'1  and  0.0042  s'1  was  used  to  generate  the 
synthetic  data.  Similar  to  what  the  model-resolution  values  indicate,  the  checkerboard-test 
result  shows  that  the  model  is  well  resolved  for  an  elongated  region  bounded  on  the  south 
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by  the  Alpine-Himalayan  orogenic  belt,  for  the  North-South  seismic  belt  in  China  and  for 
Japan.  The  best-resolved  region  both  in  terms  of  the  pattern  and  in  terms  of  the  gradient 
values  is  along  the  western  portion  of  the  Tethys  collision  belt  in  southern  Europe,  where 
the  path  density  is  also  the  highest  (Figure  64).  Again,  regions  under  northern  Eurasia, 
southeast  China,  eastern  portion  of  India  and  the  Arabian  Peninsula  are  not  well  resolved. 


Figure  71.  Checkerboard-test  result  indicating  the  resolution  power  of  the  data  used  in  the 
inversion. 

Figure  72  compares  data  residuals  using  model  predictions  from  the  a  priori  model 
and  those  from  the  final  model.  The  distribution  of  residuals  from  the  final  model  is  closer 
to  a  normal  distribution  and  less  skewed  than  the  distribution  of  residuals  from  the  a  priori 
model. 
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Figure  72.  Histograms  of  prediction  residuals  from  both  the  a  priori  model  and  from  the 
final  model. 

Figure  73  again  shows  the  tomographic  Pn  velocity-gradient  model.  Also  plotted  are 
tectonic  plate  boundaries  (Bird,  2003).  The  velocity  gradient  is  almost  universally  positive 
for  the  whole  of  Eurasia.  Much  of  the  gradient  variations  occur  along  collisional  plate 
boundaries  under  overriding  plates.  This  indicates  that  strong  velocity-gradient  variations 
and  the  presence  of  complex  tectonic  processes  involved  in  plate  collisions  may  be 
related.  The  highest  gradient  is  observed  in  the  Caucasus  region  between  the  Black  Sea 
and  the  Caspian  Sea.  Other  regions  of  higher  gradient  include  the  Scandinavia  region, 
western  and  southern  Iran,  the  Hindu  Kush  and  western  part  of  the  Himalayas.  Except  the 
Scandinavia  region,  all  other  regions  of  high  velocity  gradient  are  tectonically  active 
regions. 
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Figure  73.  The  tomographic  P„  velocity-gradient  model  and  tectonic-plate  boundaries. 

Figure  74  plots  the  RSTT  Pn  velocity-gradient  model  for  a  comparison.  As  is  indicated  by 
the  Monte  Carlo  simulation  discussed  above,  velocity  gradients  in  the  RSTT  model  are 
generally  lower.  The  mean  gradient  of  RSTT  is  0.0012  s'1  whereas  the  mean  of  the  model 
from  this  study  is  0.0044  s'1.  The  maximum  gradient  of  RSTT  is  0.0189  s'1  whereas  that  of 
our  model  is  0.0739  s'1.  There  are  also  much  less  variations  in  RSTT  than  in  our  model. 
The  difference  between  the  maximum  and  minimum  velocity  gradients  in  RSTT  is  0.0198 
s'1.  The  difference  is  0.0779  s'1  in  our  model.  The  highest  gradients  in  RSTT  concentrate 
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in  subduction  zones  in  the  western  Pacific.  There  are  also  significantly  more  negative- 
gradient  areas  in  RSTT. 


Figure  74.  RSTT  Pn  velocity-gradient  model. 


79 

Approved  for  public  release;  distribution  is  unlimited. 


5.0  SUMMARY  OF  FINDINGS 


New  2D  computations  of  Pn  geometric  spreading  have  established  that  the  amplitude 
variation  is  dependent  on  direction  of  the  path  and  that  complex  frequency  dependence  is 
expected  to  have  mixing  of  the  near-source  and  near-receiver  structures.  This  further 
weakens  the  rationale  for  using  average  power-law  representations  of  the  Pn  (and  Sn) 
geometric  spreading.  As  the  only  structure  that  will  result  in  a  power  law  behavior  for 
sure  are  those  with  negative  critical  gradients  in  the  mantle  lid  (a  possible,  but  probably 
very  rare  situation),  we  further  support  the  contention  of  Yang  (2011)  that  a  simple 
velocity  model  based  correction,  or  an  empirically  grounded  correction  is  a  more  sensible 
choice  than  a  frequency-independent  power  law. 

We  explored  2D  models  with  fine-scale  heterogeneity  to  evaluate  whether  the 
tendency  to  smooth  out  the  distance  behavior  for  1  Hz  signals  as  found  by  Avants  et  al. 
(2011)  holds  for  broadband  Pn  signals.  All  calculations,  for  suites  of  heterogeneity 
realizations,  support  some  modification  of  the  amplitude-distance  trends  for  Pn,  but  the 
frequency-dependence  of  models  with  radial  velocity  gradients  is  enhanced  by  the 
presence  of  small-scale  lithospheric  heterogeneity.  This  is  the  result  of  scattering  of  Pn 
energy  out  of  the  lid  waveguide  into  the  crust.  While  the  general  trends  of  the  amplitude 
spreading  terms  do  approach  the  power  law  decay  form  commonly  assumed,  there  is 
associated  strong  frequency  dependence,  which  is  not  commonly  assumed. 

Effects  of  crustal  thinning  and  Moho  dip  on  Pn  from  North  Korean  nuclear  tests  to 
stations  in  Japan  account  for  relatively  enhanced  Pn  amplitudes  near  1000  km  observed 
for  relatively  low  frequencies  in  Japan  compared  to  continental  structure  paths  to  stations 
in  China  near  the  same  distance.  The  lateral  transitions  from  continent  to  ocean  and  ocean 
to  island  arc  favor  transfer  of  P  energy  into  the  Pn  wavefield  near  the  sources  and  into  the 
crust  near  the  receivers.  Above  7  Hz,  the  observations  in  Japan  show  rapid  amplitude 
decrease,  which  we  attribute  to  intrinsic  attenuation  effects. 

Application  of  regionalized  Pn  spreading  corrections  based  on  RSTT  models  (derived 
exclusively  from  regional  travel  time  observations  and  velocity  models)  or  regionalized 
fitting  of  attenuation  models  to  data  provide  fairly  similar  Q(f)  estimates,  much  lower  than 
those  from  a  constant  power  law  spreading  coefficient.  While  the  choice  of  optimal 
spreading  is  not  unambiguous  in  every  case  based  on  some  measure  of  variance  reduction 
in  the  model  performance,  the  more  reasonable  attenuation  parameters  obtained  using  the 
RSTT  or  empirical  spreading  corrections  suggest  they  are  more  appropriate  than  standard 
power  law  assumptions.  This  means  that  it  is  possible  in  general  to  develop  self- 
consistent  travel  time  and  amplitude  models  for  Pn.  The  major  challenge  is  to  attain 
sufficient  data  coverage  to  actually  resolve  the  radial  gradients  in  upper  mantle  structure 
that  have  strong  effect  on  the  Pn  spreading.  When  this  is  achieved,  the  regional  model 
should  be  used  for  both  event  location  and  regional  amplitude  analysis  (and  attenuation 
model  determination),  as  the  modeling  conducted  here  demonstrates  the  reduction  of  bias 
in  the  attenuation  models  that  is  likely  to  be  achieved. 
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We  developd  a  2D  Pn  velocity-gradient  model  for  Eurasia  based  on  observed  Pn  travel 
times.  These  travel  times  were  corrected  for  crustal  and  head-wave  components  to  isolate 
the  velocity-gradient  contribution  for  path-specific  velocity-gradient  calculation.  The 
assumption  that  these  path-specific  velocity  gradients  were  the  mean  of  the  laterally 
varying  velocity  gradients  along  the  paths  was  validated  through  a  Monte  Carlo 
simulation.  A  tomographic  Pn  velocity-gradient  model  was  then  derived  using  a  linear 
Bayesian  inversion  approach. 

The  velocity-gradient  model  is  well  resolved  in  tectonically  active  regions.  The 
gradient  is  positive  almost  everywhere.  Large  gradient  variations  are  observed  mainly 
along  collisional  plate  boundaries  under  overriding  plates.  It  may  imply  a  relationship 
between  the  presence  of  active  tectonic  processes  involved  in  plate  collisions  and  the  Pn 
velocity-gradient  variation  in  the  uppermost-mantle.  This  model  can  be  used  to  predict 
path-specific  Pn  spreading  along  with  use  for  travel  time  calculation  and  event  location. 
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Highlights 

A  data-constrained,  frequency-dependent,  log-quadratic  Pn  spreading  model. 

Apparent  Pn-wave  attenuation  in  Northeast  China  and  the  Korean  Peninsula. 

Uppermost  mantle  P-wave  attenuation  model  in  the  region. 

Abstract 

We  investigate  the  geometric  spreading  and  attenuation  of  seismic  Pn  waves  in 
Northeast  China  and  the  Korean  Peninsula.  A  high-quality  broadband  Pn-wave  dataset 
generated  by  North  Korean  nuclear  tests  is  used  to  constrain  the  parameters  of  a 
frequency-dependent  log-quadratic  geometric  spreading  function  and  a  power-law  Pn  Q 
model.  The  geometric  spreading  function  and  apparent  Pn  wave  Q  are  obtained  for 
Northeast  China  and  the  Korean  Peninsula  between  2.0  and  10.0  Hz.  Using  the  two-station 
amplitude  ratios  of  the  Pn  spectra  and  correcting  them  with  the  known  spreading  function, 
we  strip  the  contributions  of  the  source  and  crust  from  the  apparent  Pn  Q  and  retrieve  the 
P-wave  attenuation  information  along  the  pure  upper  mantle  path.  We  then  use  both  Pn 
amplitudes  and  amplitude  ratios  in  a  tomographic  approach  to  obtain  the  upper  mantle  P- 
wave  attenuation  in  the  studied  area.  The  Pn-wave  spectra  observed  in  China  are 
compared  with  those  recorded  in  Japan,  and  the  result  reveals  that  the  high-frequency  Pn 
signal  across  the  oceanic  path  attenuated  faster  compared  with  those  through  the 
continental  path. 

Key  words 

Pn;  Geometric  spreading;  Attenuation;  Uppermost  mantle  P-wave  Q;  North  Korea 
nuclear  tests 
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1  Introduction 


The  seismic  Pn  wave  typically  appears  as  the  first  arrival  in  regional  seismograms  at 
distances  between  200  and  2000  km.  Unlike  other  regional  phases,  such  as  Pg  and  Lg 
waves,  the  Pn  wave  is  typically  less  affected  by  the  crustal  structures  or  contaminated  by 
prior  phases  and  is  therefore  widely  applied  for  seismic  location,  magnitude  and  yield 
estimation  as  well  as  for  investigating  seismic  source  process.  The  Pn-wave  attenuation 
measurement  is  an  important  issue  as  well  not  only  because  obtaining  accurate  source 
parameters  requires  correction  for  the  amplitude  decay  but  also  attenuation  itself  can  be 
useful  for  characterizing  the  upper  mantle  properties.  Pn  and  Sn  waves  cannot  be  simply 
explained  as  headwaves  developed  along  a  flat  Moho  discontinuity.  Their  actual 
propagation  processes,  particularly  their  amplitude  variations  in  a  spherical  earth,  are 
much  more  complicated.  Cerveny  and  Ravindra  [1971]  and  Hill  [1973]  theoretically 
investigated  the  behavior  of  Pn  waves  in  a  spherical  Earth  model  and  regarded  this  phase 
as  an  interference  of  multiple-diving  waves  reflected  from  the  underside  of  the  Moho 
discontinuity.  Sereno  and  Given  [1990]  studied  Pn  waves  in  flat  and  spherical  earth 
models  and  found  that  Earth’s  sphericity  alone  causes  a  significant  departure  in  Pn 
geometrical  spreading  and  that  the  phenomenon  is  frequency  dependent.  Based  on 
numerical  simulations,  Yang  et  al.  [2007]  proposed  a  log-quadratic  model  for  Pn  and  Sn 
geometric  spreading  functions  to  accommodate  the  Earth’s  sphericity.  Compared  with  the 
traditional  power-law  model,  their  spreading  model  included  9  parameters  to  address  both 
distance-  and  frequency-dependencies. 

Many  attempts  have  been  targeted  at  using  observed  data  to  constrain  Pn-wave 
attenuation  and  geometrical  spreading  [Chun  et  al.,  1989;  Sereno  et  al.,  1988;  Xie,  2007; 
Xie  and  Patton,  1999;  Zhu  et  al.,  1991].  The  underlying  difficulties  are  that  attenuation 
and  geometrical  spreading  are  tightly  coupled  and  both  are  dependent  on  distance  and 
frequency.  In  addition,  the  observed  Pn  amplitudes  are  highly  scattered  due  to  their 
narrow  sampling  of  radiation  patterns,  inaccurate  source  locations,  and  pronounced 
sensitivity  to  the  uppermost  mantle  structures.  Different  strategies  were  adopted  to 

mitigate  this  difficulty.  By  assuming  a  frequency- independent  spreading  function  of  A  13 , 
where  A  is  the  distance,  Sereno  et  al.  [1988]  investigated  the  Pn-wave  Q  in  Scandinavia 
and  obtained  a  Qo  (1  Hz  Q)  of  325.  Using  the  same  spreading  function,  Xie  and  Patton 
[1999]  obtained  a  Qo  of  364  for  central  Asia,  and  Xie  [2007]  observed  low  Pn  Q  under  the 
north-central  Tibetan  plateau  between  0.3  and  10.0  Hz.  Chun  et  al.  [1989]  determined  the 
high-frequency  Pn  attenuation  in  eastern  Canada  based  on  a  frequency-dependent  Pn 

spreading  z\"(2-2+0-02'/)  between  3.0  and  15.0  Hz.  Zhu  et  al.  [1991]  simultaneously  estimated 
both  frequency-dependent  geometric  spreading  and  Q  for  Pn  waves  in  the  Canadian 
shield. 

Others  attempted  to  estimate  the  Pn-  or  Sn-wave  propagation  efficiency  without 
apparently  addressing  the  geometric  spreading.  Calvert  et  al.  [2000]  measured  the 
propagation  efficiency  by  separating  the  attenuation  into  the  average  portion  and  the 
perturbation  portion,  avoiding  defining  the  absolute  attenuation  measurement.  Using  the 
ratio  of  the  Sn  amplitude  to  the  Pg  coda  amplitude,  Barron  and  Priestley  [2009]  presented 
the  frequency-dependent  propagation  efficiency  of  the  Sn  wave  over  the  Tibetan  Plateau. 
By  stacking  the  densely  distributed  US  Array  data,  Buehler  and  Shearer  [2013]  calculated 
the  Sn  propagation  efficiency  and  identified  highly  attenuating  regions  in  the  Western 
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United  States.  Another  strategy  involves  building  a  geometric  spreading  model  followed 
by  separating  the  attenuation  from  the  spreading.  Yang  et  al.  [2007]  proposed  a  log- 
quadratic  geometric  spreading  function  based  on  numerical  simulations  in  a  radially 
symmetric  earth  model.  Avants  et  al.  [2011]  further  discussed  the  contribution  of 
scattering  to  this  spreading  model.  Applying  this  log-quadratic  spreading  model  to  Eurasia 
data,  Yang  [2011]  obtained  a  more  reasonable  Pn  Q  between  distances  of  200  and  1000 
km  and  frequencies  of  0.5  and  10.0  Hz. 

Between  2006  and  2013,  North  Korea  conducted  three  underground  nuclear  tests. 
Compared  with  natural  earthquakes,  these  explosive  sources  have  accurate  epicenter 
locations  and  source  depths,  simple  and  approximately  identical  source  time  functions, 
and  most  importantly  virtually  isotropic  radiation  patterns.  Seismic  networks  recorded 
abundant  regional  phases  from  these  events,  including  the  Pn  wave,  in  Northeast  China, 
the  Korean  Peninsula  and  Japan,  either  across  the  continental  or  oceanic  paths  [e.g.,  Chun 
et  al.,  2009;  Hong  and  Rhie,  2009;  Murphy  et  al.,  2013;  Richards  and  Kim,  2007;  Wen  and 
Long,  2010;  Zhao  et  al.,  2008,  2012,  2014], 

In  this  study,  we  use  this  high-quality  dataset  to  investigate  Pn-wave  propagation  in 
Northeast  China  and  the  Korean  Peninsula,  including  its  geometrical  spreading  and  the 
apparent  Pn  attenuation.  An  inversion  method  is  used  to  formally  separate  the  spreading 
and  attenuation  using  the  observed  spectral  amplitudes.  Frequency  dependent  parameters 
for  a  log-quadratic  spreading  model  in  Northeast  China  and  the  Korean  Peninsula  are 
derived  from  the  data.  Using  two-station  amplitude  ratios  and  the  known  spreading  function, 
we  strip  the  effect  of  the  crust  leg  from  the  apparent  Pn-wave  attenuation  to  obtain  the  P- 
wave  attenuation  along  the  pure  upper  mantle  path.  Combining  both  single-station  amplitudes 
and  two-station  amplitude  ratios,  an  uppermost  mantle  P-wave  attenuation  tomography 
approach  is  proposed.  Using  Pn  observations  across  Northeast  China  and  those  through 
the  Japan  Sea,  the  Pn  spectral  amplitudes  across  the  continental  and  oceanic  paths  are 
compared. 

2  Data 

2.1  Regional  dataset 

On  9  October  2006,  25  May  2009,  and  12  February  2013,  North  Korea  conducted 
three  successive  nuclear  tests  near  the  China-Korea  border,  and  their  body  wave 
magnitudes  were  4.2,  4.7  and  5.1,  respectively,  as  reported  by  the  United  States 
Geological  Survey  (USGS).  Hereafter,  these  events  are  referred  to  as  NKT1,  NKT2  and 
NKT3,  and  their  parameters  are  provided  in  Table  1.  Broadband  digital  seismograms  from 
these  nuclear  tests  are  collected.  Figure  1  depicts  the  locations  of  the  North  Korean 
nuclear  test  site  (NKTS)  and  the  297  broadband  stations  used  in  this  study.  Among  them, 
188  permanent  stations  are  from  the  national  and  provincial  networks  under  the  China 
Earthquake  Administration  (CEA)  and  the  China  National  Digital  Seismic  Network 
(CNDSN)  operated  by  the  China  Earthquake  Networks  Center  (CENC)  since  December 
2000.  Seven  permanent  stations  are  from  the  Global  Seismographic  Network  (GSN), 
which  has  been  operated  by  the  USGS  and  the  Incorporated  Research  Institutions  for 
Seismology  (IRIS)  consortium  from  May  1994  to  present.  Two  portable  seismic  arrays, 
including  the  6  station  SASK  in  South  Korea  and  the  35  station  NECsaids  in  Northeast 
China,  originally  for  investigating  deep  structures,  recorded  NKT2  and  NKT3  with  high 
signal-to-noise  ratios.  The  remaining  61  stations  are  from  F-NET  in  Japan.  These  stations 
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are  equipped  with  broadband  instruments  with  nearly  flat  velocity  responses  between  0.03 
and  8.0  Hz  or  wider,  and  their  sampling  rates  vary  among  20,  40,  50,  and  100  points  per 
second.  The  distances  from  the  NKTS  to  these  stations  are  between  143  and  1930  km, 
within  which  the  Pn  phase  is  well-developed. 

We  plot  regional  waveforms  in  a  record  section  to  investigate  the  Pn-wave  group 
velocity  in  this  area.  Instead  of  plotting  waveforms  themselves,  we  plot  normalized  waveform 
energy.  Figure  2a  presents  a  vertical-component  velocity  seismogram  from  NKT3  recorded 
by  station  LN.LYA  at  a  distance  of  813.0  km,  where  a  bandpass  filter  between  0.1  and 
10.0  Hz  is  applied.  Figure  2b  shows  the  normalized  waveform  energy.  Figure  2c  presents 
the  waveform  energy  versus  distance,  in  which  several  regional  phases  can  be  identified. 
The  first  arrival  is  the  Pn  wave,  which  has  a  group  velocity  of  approximately  8.1  km/s. 
The  Pg-  and  Lg-waves  can  be  traced  at  group  velocities  of  6.0  and  3.5  km/s,  respectively. 

2.2  Pn  amplitude  measurements 

The  Pn  wave  is  typically  the  first  arrival  at  regional  distances,  and  its  amplitude  can 
be  measured  within  a  time  or  group  velocity  window.  For  example,  a  fixed  5-second  time 
window  was  used  in  Scandinavia  by  Sereno  et  al.  [1988],  a  4.3-second  window  was  used 
in  Eastern  Canada  by  Zhu  et  al.  [1991],  and  a  4.5-second  window  was  used  in  Tibet  by  Xie 
[2007].  Considering  the  slightly  dispersive  properties  of  the  Pn  wave,  some  authors  have 
used  a  group  velocity  window  between  8.2  and  7.6  km/s  to  measure  the  amplitude  [e.g., 
Al-Damegh  et  al.,  2004;  Mcnamara  et  al.,  1997;  Reese  et  al.,  1999].  In  this  study,  we  use 
the  vertical-component  seismogram  and  a  0.7  km/s  group-velocity  window  around  the 
IASP91  arrival  time  to  measure  the  Pn  amplitude.  Figure  3  briefly  illustrates  the  data 
processing  process.  The  vertical-component  seismogram  recorded  at  station  LN.LYA 
from  NKT3  is  presented  in  Figure  3a,  where  the  gray-shaded  areas  indicate  the  0.7  km/s 
group  velocity  window  for  measuring  Pn  and  the  pre-P-arrival  window  for  sampling  the 
noise.  The  windowed  signal  is  plotted  in  Figure  3b,  where  10%  cosine  tapers  are  applied 
at  both  ends.  Following  Zhao  et  al.  [2010,  2013b],  we  select  the  noise  series  in  an  equal- 
length  time  window  before  the  first  arrival.  Figure  3c  presents  the  calculated  Fourier 
spectra  of  the  Pn  and  noise  series  between  0.3  and  15.0  Hz.  The  Pn  spectral  amplitudes  are 
obtained  at  44  frequencies  distributed  log-evenly  between  0.3  and  15.0  Hz  and  corrected 
for  the  noise  (Figures  3d  and  3e).  (for  details,  see  Zhao  et  al.,  2013a;  Zhao  et  al.,  2010). 


3  Pn-wave  geometric  spreading  and  attenuation 
3.1  Modeling  of  the  Pn  spectrum 

The  Pn-wave  spectrum  can  be  expressed  as  [Sereno  et  al.,  1988;  Xie,  2007;  Xie  and 
Patton,  1999]: 

^(/)  =  S(/)-G(A,/)-r(A,/).  />(/).  r(f),  (1) 

where  /  is  the  frequency,  -d(f)  is  the  observed  spectral  amplitude,  P{f)  is  the  site 

response,  and  r{f  )  is  the  random  amplitude  effect.  S(/)  is  the  Pn  source  spectrum, 

which  is  given  by  the  following  simplified  explosion  source  function  [Hong,  2013; 
Mueller  and  Murphy,  1971;  Sereno  et  al.,  1988;  Xie  and  Patton,  1999]: 


S(f)  =  S0 


1+(1_2/?)T+/?!T' 


-1/2 


(2) 
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where  S0  is  the  long  period  source  spectral  level,  fc  is  the  comer  frequency,  and  J3 
controls  the  amount  of  overshot.  For  a  Poisson  medium,  /]  =  0.75.  For  a  near-surface 
explosion,  S0  =  M 0l / Anpor1 ,  where  M0  is  the  seismic  moment,  and  p  and  a  are  the 

density  and  P-wave  velocity  in  the  source  region  [Sereno  et  al.,  1988;  Stevens  and  Day, 
1985],  which  in  this  study  are  set  to  2.7  g/cm3  and  5.5  km/s  [Hong  et  al.,  2008;  Jih,  1998], 

respectively.  In  Equation  (1),  the  geometric  spreading  factor  G(A,/)  is  a  function  of  the 
epicenter  distance  A  and  frequency  / ,  and  its  details  are  provided  in  the  next  section.  In 
Equation  (1), 


r(A,/)  =  exp 


(3) 


is  the  attenuation  factor,  where  V  is  the  Pn-wave  group  velocity,  and  will  be  treated  as  a 
regional  constant, 


ds 


’Q(f) 


(4) 


is  the  integral  of  attenuation  over  the  great  circle  path.  In  addition,  Q  (/)  is  the  apparent 

Pn-wave  Q  and  assumed  can  be  expressed  by  a  power-law  model  [e.g.,  Sereno  et  al., 
1988;  Xie,  2007]: 


Q{f)  =  QJ\ 


(5) 


where  Q0  and  r/  are  the  apparent  Pn-wave  Q  at  1  Hz  and  its  frequency-dependence, 
respectively. 


3.2  Pn-wave  geometric  spreading  function 

The  power-law  spreading  function  linearly  relates  the  logarithmic  amplitude  with  the 
logarithmic  distance.  To  account  for  the  complex  spreading  relation  obtained  from  the  Pn- 
wave  simulation  in  a  spherical  Earth,  Yang  et  al.  [2007]  extended  the  power-law  relation 
by  adding  a  second-order  term.  Considering  that  Pn  spreading  is  frequency  dependent, 
these  researchers  also  included  a  quadratic  form  for  its  frequency  dependency.  We  will 
call  the  model  the  log-quadratic  spreading  model,  which  exhibits  the  following  form 
[Yang  et  al.,  2007]: 

G(A,/)  =  10«/)a4(/)Ios“4-«/)j  (6) 

where  coefficients  %n{f)  are  dependent  on  the  logarithm  of  the  frequency 


4  (/)  =  4i  logfo  (/)  +  4,2  lo§io  (/)  +  4,3  >  ( n  =1 ,  2,  3).  (7) 

Based  on  a  spherical  earth  model  composed  of  a  40-km  thick  crust  and  an  underlying 
constant  upper  mantle,  Yang  et  al.  [2007]  used  numerical  simulation  to  determine  the 
coefficients  as 


-0.217  1.79  3.16 
-1.94  8.43  18.6 

-3.39  9.94  20.7 


(8) 


Equation  (6)  is  likely  a  more  flexible  model  for  the  Pn  geometric  spreading  than  the 
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simple  power-law  model,  but  its  coefficients  can  vary  for  regions  with  different  velocities, 
crust  thickness,  and  upper  mantle  lid  structures.  Therefore,  we  adopt  equation  (6)  for 
Northeast  China  and  the  Korean  Peninsula  but  fine-tune  its  parameters  using  observed  Pn 
data. 

3.2.1  Prior  information  and  model  space 

Table  2  summarizes  the  model  parameters  to  be  searched  for,  including  four 
parameters  for  sources,  nine  geometric  spreading  coefficients,  and  two  parameters  for  the 
apparent  Pn-wave  Q.  We  estimate  the  variation  range  of  these  parameters.  For  the  three 
North  Korean  nuclear  tests,  their  seismic  moment  M0  and  the  comer  frequency  fc  can  be 

calculated  from  their  body-wave  magnitude  mb  using  the  empirical  relationship 
established  by  Taylor  et  al.  [2002],  We  set  ±10%  perturbations  for  mb  to  obtain  the 
variation  ranges  for  M0.  The  source  comer  frequency  fc  can  be  calculated  from  the 
seismic  moment  using  the  empirical  relation  log M()  =  a +  b- log  fc  from  Xie  and  Patton 
[1999],  From  the  Pn  data  recorded  for  explosions  in  central  Asia,  these  researchers  found 
a  =  18.34  and  b  =  4.73 .  There  is  likely  a  tradeoff  between  M0  and  a ,  so  we  fix  a  =  18.34 

in  our  model.  The  slope  b  also  exhibits  certain  regional  dependence,  so  we  allow  ±30% 
variation  for  this  parameter.  The  log-quadratic  Pn  spreading  model  derived  by  Yang  et  al. 
[2007]  is  expected  to  be  transferable  to  other  regions  with  proper  modification  of  its 
parameters.  Thus,  we  impose  ±10%  perturbation  for  all  9  coefficients.  Sereno  et  al.  [1988] 

obtained  apparent  attenuation  of  Q(f)  -  325/0'48  in  Scandinavia  using  1  to  15  Hz  Pn 
spectra.  Sereno  [1990]  investigated  Pn  amplitude  spectra  in  eastern  Kazakhstan  in  central 
Asia  and  obtained  a  similar  power-law  Q  model  of  Q{f  )  =  300/° 5 .  Chun  et  al.  [1989] 

suggested  that  a  frequency-dependent  Q  model,  Q(f)  -  255/0'74 ,  can  explain  the 
observed  Pn  spectra  in  the  Canadian  shield  between  3  and  10  Hz.  Xie  [2007]  used  0.6  to 
5.0  Hz  Pn  spectra  to  obtain  an  average  Q  model,  Q{f)  =  278 /014 ,  for  the  Tibetan  plateau. 
Yang  [201 1]  investigated  Pn  geometric  spreading  based  on  observations  in  Asia,  wherein 
he  observed  that  the  Pn  Q  approximately  follows  a  power-law  function  Q(f)  =  204 /° 58 
for  frequencies  greater  than  2  Hz.  Using  these  results  as  prior  information,  we  limit  the 
range  of  Q0  and  //  from  200  to  600  and  from  -0.2  to  1.2,  respectively.  All  prior 
information  and  model  parameter  spaces  are  summarized  in  Table  2. 

3.2.2  Pn  geometric  spreading  model  in  Northeast  China  and  the  Korean  Peninsula 

Because  of  the  potentially  different  Pn  propagation  characteristics  of  continental  and  oceanic 
paths,  and  because  the  continental  dataset  covers  a  wide  and  relatively  even  distance  range,  we  use 
only  continental-path  Pn  amplitude  measurements  in  the  following  parameter  searching  and 
tomography  analysis.  We  also  restrict  our  analysis  to  amplitudes  within  2.0  to  10.0  Hz  due  to 
apparently  different  Pn  attenuation  behavior  (see  Figure  7d)  beyond  this  frequency  range. 
Simulated  annealing  [Kirkpatrick  et  al.,  1983],  a  non-exhaustive  global  optimization 
algorithm,  is  used  to  estimate  the  parameters  in  model  space.  This  method  has  been  widely 
applied  in  geophysical  modeling  [e.g.,  Iritani  et  al.,  2014;  Kirkpatrick,  1984;  Zhao  et  al., 
1996].  We  perform  the  parameter  searching  by  minimizing  the  L2  norm  misfit  function  of  the  difference 
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between  observed  and  synthetic  Pn  amplitudes.  Figure  4  illustrates  the  Pn  spectral  amplitudes  at 
selected  frequencies,  after  correcting  for  the  source  excitation  functions.  Different 
symbols  indicate  that  the  samples  are  from  NKT1  (triangles),  NKT2  (circles)  and  NKT3 
(crosses).  The  solid  blue  lines  are  the  best-fit  Pn  amplitudes  from  data  between  2.0  and 
10.0  Hz,  and  the  shaded  areas  are  their  standard  deviations  of  the  fit  calculated  using  the 
bootstrap  method  [Efron,  1983],  The  frequency  f  number  of  samples  N,  apparent  regional 
average  Pn  Q,  and  correlation  coefficient  R  are  labeled  in  each  panel. 

The  inverted  Pn  models  are  listed  in  Table  2  and  illustrated  in  Figure  5.  Figures  5a 
and  5b  presents  the  Pn  geometric  spreading  curves  versus  distance  at  2.1  and  9.1  Hz, 
respectively.  Figures  5c  and  5d  present  these  curves  versus  frequency  at  distances  of  500 
and  1,200  km,  respectively.  The  solid  lines  indicate  the  results  obtained  here,  and  the 
dashed  lines  represent  the  model  proposed  by  Yang  et  al.  [2007],  The  two  spreading 
models  exhibit  similar  overall  shapes  that  first  decrease  and  then  increase  with  increasing 
distance.  The  minimums  appear  at  shorter  distances  for  higher  frequencies.  The  models 
are  also  positively  frequency  dependent,  at  least  within  the  investigated  frequency  band. 
Figure  5e  illustrates  the  inverted  source  excitation  functions,  where  the  shaded  areas  are 
the  standard  deviations.  Figure  5f  presents  the  average  regional  Pn  Q  versus  the  frequency. 

4  Upper  mantle  P-wave  attenuation 

Pn-wave  attenuation  is  an  apparent  value  because  the  Pn  Q  observed  along  the  great 
circle  is  mixed  with  those  traveling  in  the  uppermost  mantle  and  crust,  with  the  latter 
propagating  with  a  nonzero  dipping  angle.  With  the  known  Pn-wave  geometric  spreading 
function,  we  can  strip  the  effect  of  crust  legs  and  obtain  the  upper  mantle  P-wave 
attenuation. 


4.1  Single-station  data 

Figure  6  presents  the  Pn  propagation  path  between  the  NKTS  and  station  LN.LYA, 
calculated  using  the  CRUST  1.0  model  [ Laske  et  al.,  2013],  The  refraction  points  E  and 
F  separate  the  ray  path  into  three  sections,  including  the  crust  leg  at  the  source  side,  the 

uppermost  mantle  leg,  and  the  crust  leg  at  the  station  side.  From  equation  (4),  B (A,/) 
can  be  expressed  as 

.  r  ds  r  ds  r  ds 

/)  =  — +  - r  +  — »  (9) 

'  ]aeQs  EF  Q(x,y,f)  }™Qr 

where  Qs  and  Qr  are  the  source  and  station  side  crust  P-wave  Q,  and  Q(x,  yj )  is  the 


uppermost  mantle  P-wave  Q,  which  is  a  function  of  the  frequency  and  location  (v,y) 

between  E  and  F  .  By  substituting  Equation  (9)  into  (3),  the  Pn  attenuation  factor  can  be 
expressed  as 


F(A,/)  =  Y  (AEJ)  Y(EF,f )  •  r,  (FB,f) 


r  ds 

ds 

nf 

r  ds 

«u  ' 

' AE~QS 

a2  >EFQ(x,y,f ) 

<*1  r  ’ 

>FBQ_ 

(10) 


where  als  and  alr  are  the  source  and  station  side  P-wave  velocity,  a2  is  the  upper  mantle 
P-wave  velocity.  From  equations  (1)  and  (10),  we  have 
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A (f)G~‘  (A, /) ■  S-'  (/)  =  [r,  (AE,f)T,  (FB,f)\  ■  exp 


Of, 


f  - 


rA 


,(H) 


Q{x,y,f ) 

where  G(A,/)  and>S(/)  are  known  from  Section  3,  and  A(f)G  '(A,/)^  ’(/)  is  the 

spreading  and  source  corrected  single- station  data.  T s{AE,f}-Tr(FB,f}  is  the 

attenuation  from  the  crust  legs.  We  assume  it  is  an  unknown  constant  over  the  studied 
region  and  can  be  inverted  for  from  the  data.  It  is  expected  that  the  upper  mantle  leg 
dominates  the  Pn  propagation,  and  the  error  in  the  crust  legs  can  be  neglected.  We  will 

also  neglect  factors  P(f)  and  r  (/)  caused  by  the  site  effect  and  the  random  fluctuations. 


4.2  Two-station  data  and  Pn  amplitude  ratio 

Assuming  that  the  two  stations  are  located  at  different  epicenter  distances  but  within 
the  similar  azimuth  direction  from  the  NKTS,  we  calculate  their  amplitude  ratio  as  follows 
[e.g.,  Xie  et  al.,  2004]: 

A  MU  G(Mf)  EiEull r-  (FB-f)  1  (lln (/) 

*  4 (/)  G(\,f)r,(AE„f)r(EFlJ)rr(FB„f)  P,(f)-r,(f)  ;  (12) 

where  Aj  (/) and  A.  (/)  are  the  Pn  amplitude  spectra  observed  at  stations  i  and  j ,  the 
source  terms  have  been  canceled,  and  G(Ay, /)/(?( A.,/)  is  the  ratio  between  known 

geometrical  spreading  functions  at  these  distances.  TsyAEt,f)  and  Tv (dfT./j  represent 
the  source  side  crust  legs.  For  observations  from  the  same  source,  their  ratio  should  be 
equal  to  unity.  T;  (FBnf)  and  Tr  represent  the  station  side  crust  legs.  As  above 

mentioned,  we  assume  these  values  are  approximately  the  same,  and  their  ratio  approaches 
unity.  We  also  neglect  the  ratio  Pj  (/)  •  Tj  (f)/Pj  (/)  •  (/) .  After  these  treatments  and  by 
introducing  equation  (10),  equation  (12)  becomes 


A„(f) 


G(A„/) 

C(Aj,f) 


exp 


a , 


ds 


Q(x,y,f) 


(13) 


c  F 

where  the  integral  \j Q(x,y,  f)ds  is  from  the  refraction  point  Fj  to  F j  over  the  upper 

mantle  path.  The  left  hand  side  of  equation  (13),  the  amplitude  ratio  corrected  by  the 
known  geometrical  spreading  function,  is  directly  linked  to  the  accumulated  attenuation 
over  the  pure  upper  mantle  path.  Thus,  the  two-station  amplitude  ratio  strips  the  source 
term  and  contributions  from  the  crust  legs  near  the  source  and  station. 

Figures  7a  to  7c  illustrate  the  Pn  amplitude  ratios  after  correction  for  the  spreading 
function  versus  the  interstation  distance  at  3.0,  5.2  and  9.1  Hz.  The  solid  lines  are  linear 
regressions,  and  their  slopes  provide  the  regional  average  upper  mantle  P-wave  Q.  The 
circles  in  Figure  7d  represent  upper  mantle  P-wave  Q  versus  frequency.  They  show  a 
nearly  linear  relation  between  2.0  and  10.0  Hz  (shaded  area),  within  which  a  power-law  Q 

model,  QP  =17 6 /0  48 ,  for  the  upper  mantle  P  wave  can  be  obtained.  However,  outside  this 
frequency  band,  the  Q  versus  frequency  behaves  quite  differently.  In  our  Pn  Q  model,  we 
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assume  a  power-law  frequency  dependence.  The  results  demonstrate  that  this  assumption 
may  only  be  valid  within  the  2.0  to  10.0  Hz  band  in  the  studied  region. 

4.3  Upper  mantle  P-wave  Q  tomography 

An  inversion  system  for  the  upper  mantle  P-wave  Q  can  be  created  either  from  the 
single-station  data,  the  two-station  data,  or  by  combining  both.  One  might  think  that  the 
two  datasets  are  from  the  same  group  of  observations  and  should  contain  exactly  the  same 
information,  but  certain  differences  still  exist.  The  advantage  of  using  the  two-station  data 
is  that  the  amplitude  ratio  eliminates  the  source  effect  at  the  data  processing  stage, 
reducing  the  tradeoff  between  the  source  term  and  attenuation  in  the  inversion.  However, 
the  slight  azimuth  difference  between  two  stations  may  result  in  additional  errors.  Most 
importantly,  the  two-station  data  must  match  the  source-station  geometry,  which  often 
results  in  less  useful  rays  and  shorter  ray  length  compared  with  the  single- station  method. 
This  effect  is  particularly  severe  for  the  heavily  attenuated  region  and  for  high 
frequencies,  where  near-station  information  is  vital.  Thus,  combining  the  single-  and  two- 
station  data  improves  data  coverage  and  avoids  serious  tradeoffs. 

Similar  to  the  Lg-wave  Q  tomography  [e.g.,  Zhao  et  al.,  2010,  2013b],  we  linearize 
equations  (11)  and  (13)  by  taking  the  logarithm  to  combine  both  the  single-  and  two-station 
data  for  inversion.  In  our  tomography,  the  initial  Q  model  is  a  constant  model  resulted  from  two- 
station  data  analysis,  as  shown  in  Figure  7.  The  constant  Q  model  is  also  used  to  create 
checkerboard-shaped  model  by  superimposing  positive  and  negative  perturbations.  A  broadband 
uppermost  mantle  P-wave  attenuation  model  in  Northeast  China  and  the  Korean  Peninsula 
at  18  discrete  frequencies  between  2.0  and  10.0  Hz  is  obtained.  As  examples,  Figure  8a 
presents  the  P-wave  Q  distribution  at  3.0  Hz,  and  Figure  8b  presents  the  broadband  Q, 
which  is  calculated  by  average  logarithmic  Q  between  2.0  and  10.0  Hz.  The  major 
geology  blocks,  including  the  Songliao  Basin  (SB),  Bohai  Bay  Basin  (BB)  and  Changbai 
Mountains  (CM),  are  illustrated  in  these  figures.  The  northern  SB,  the  areas  between  the 
southern  SB  and  BB,  and  the  southeastern  CM  are  characterized  by  low  Q  anomalies, 
whereas  areas  between  SB  and  CM  are  high  Q  regions.  A  number  of  volcanoes  in  this 
region  were  previously  investigated  using  seismic  velocity  tomography  [e.g.,  Duan  et  al., 
2009;  Lei  and  Zhao,  2005;  Zhao  et  al.,  2009;  Zhao  et  al.,  2011].  Prominent  low  P-wave 
velocity  anomalies  have  been  found  in  the  crust  and  upper  mantle  beneath  these  volcanoes 
[Zhao  et  al.,  2009].  These  anomalies  are  consistent  with  the  strong  P-wave  attenuation  in 
the  uppermost  mantle  observed  in  this  study.  Figures  8c  and  8d  illustrate  the  3.0-Hz  ray 
coverage  and  the  checkerboard  test  for  the  2.0°  x  2.0°  resolution. 

5  Comparison  of  Pn-wave  amplitudes  through  continental  and  oceanic  paths 

The  ray  paths  from  NKTS  to  Japan  cross  the  Japan  Sea,  where  the  oceanic  Moho 
depth  is  approximately  1 1  to  15  km  [Laske  et  al.,  2013].  The  continental  crust  and  oceanic 
crust  have  different  thicknesses  and  underlying  upper  mantle  P-wave  velocities,  and  the 
transition  zone  between  them  may  also  influence  Pn  attenuation  and  spreading.  In  Figure 
1,  the  data  observed  via  continental  path  surround  the  NKTS  for  an  azimuth  range  of 
approximately  160°  and  cover  epicenter  distances  from  150  to  1300  km,  forming  the 
dataset  for  investigating  the  Pn  spreading  function  and  attenuation  in  Northeast  China  and 
the  Korean  Peninsula.  The  data  crossing  the  oceanic  path  are  mostly  recorded  by  F-NET 
stations  in  Japan.  The  data  cover  an  azimuth  range  of  approximately  120°  but  are  mostly 
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recorded  at  distances  of  approximately  1,000  km,  thus  preventing  them  from  being  used 
for  an  independent  geometric  spreading  function  for  the  oceanic  path. 

Given  that  the  explosion  sources  are  virtually  isotropic  and  the  observations  are 
nearly  completely  surrounding  the  source,  the  data  still  provide  an  excellent  opportunity  to 
investigate  the  Pn  amplitudes  and  their  frequency  dependency  on  the  continental  and 
oceanic  paths.  To  avoid  ambiguity,  we  make  a  straight  forward  comparison  with  minimum 
data  processing.  Given  that  the  oceanic  data  are  mostly  recorded  at  1 ,000  km,  we  collect 
the  Pn  spectral  amplitudes  between  distances  800  and  1 ,200  km  (with  a  nominal  epicentral 
distance  of  1,000  km).  We  remove  the  Pn-wave  source  excitation  functions  and  present 
the  result  in  Figure  9  as  a  function  of  the  azimuth.  The  triangles,  circles  and  crosses  are  the 
directly  measured  spectral  amplitudes  for  NKT1,  NKT2  and  NKT3.  In  addition,  black, 
blue  and  red  colors  indicate  0.8-,  7-  and  10.0-Hz  data,  respectively.  Solid  circles  with  error 
bars  indicate  the  mean  values  and  standard  deviations  obtained  within  a  30°  azimuth 
window.  Prominent  differences  can  be  observed  for  the  Pn  waves  crossing  different  paths. 
Between  60°  and  180°,  the  Pn  signals  crossing  the  oceanic  path  are  strongly  frequency 
dependent,  with  extremely  low  high-frequency  amplitudes.  The  Pn  waves  across  the 
continental  path  can  be  divided  into  two  groups.  Between  230°  and  280°  with  the  direction 
towards  Northern  China,  the  low-frequency  contents  (0.8  and  7  Hz)  are  similar  to  that 
from  the  oceanic  path,  but  the  high-frequency  (10  Hz)  propagation  is  much  more  efficient 
than  that  from  the  oceanic  path.  Within  the  azimuth  range  of  315°  to  30°  towards 
Northeast  China,  the  Pn  spectra  are  less  dependent  on  the  frequency.  In  general,  the 
intermediate  frequency  (7  Hz)  spectra  are  less  affected  by  different  paths,  but  the  low-  and 
high-frequency  contents  exhibit  apparent  variations  between  the  continental  and  oceanic 
paths.  Because  no  data  processing,  other  than  the  noise  and  source  excitation  function 
removal  is  involved,  the  observations  are  relatively  reliable.  The  spreading  function  tends 
to  raise  the  high-frequency  signal,  whereas  attenuation  tends  to  reduce  the  high-frequency 
contents.  The  detailed  structure  of  transition  zones  may  also  affect  the  frequency 
dependence.  Additional  observational  and/or  numerical  investigations  are  required  to 
distinguish  contributions  from  individual  mechanisms. 

6  Discussion 

To  determine  the  Pn  wave  geometric  spreading  function  and  the  attenuation,  several 
assumptions  are  adopted,  e.g.,  isotropic  source  radiation,  a  log-quadratic  spreading  model, 
and  a  power-law  Pn  Q  model.  Given  that  both  geometric  spreading  and  attenuation  are 
frequency  dependent,  tradeoffs  between  them  may  exist.  Certain  clues  (e.g.,  Figure  7d) 
indicate  that  a  power-law  Q  model  may  oversimplify  the  frequency  dependency  of  the 
attenuation  in  a  broad  frequency  band.  In  Section  3,  we  obtain  the  regional  average  Pn  Q 
using  the  single-station  data,  and  we  obtain  the  regional  average  upper  mantle  P  wave  Q 
using  the  two-station  data  in  Section  4.  However,  because  the  single-station  data  extend  to 
a  larger  region  compared  with  the  two-station  data,  the  two  average  Qs  do  not  necessarily 
cover  the  same  region.  In  the  upper  mantle  P  wave  attenuation  tomography,  all  seismic 
rays  come  from  the  same  epicenter  and  there  are  no  crossing  rays.  This  forms  an 
unfavorable  geometry  for  a  high-quality  tomography.  Therefore,  the  result  is  limited  and 
primarily  used  for  introducing  the  technique. 

The  Pn-wave  has  a  typical  group  velocity  window  between  8.2  and  7.6  km/s  for 
continental  paths.  Selecting  a  window  and  the  component  to  measure  the  Pn  spectra  is  an 
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important  issue.  Sereno  et  al.  [1988]  tested  different  window  lengths  of  5  to  20  s  and 
suggested  that  the  spectra  are  insensitive  to  the  window  length.  In  fact,  their  shortest  5- 
second  window  generated  approximately  the  same  Pn  measurements  for  low-magnitude 
earthquakes  with  epicentral  distances  between  200  and  1400  km.  To  investigate  the 
stability  of  the  Pn-wave  amplitude  measurement,  we  tested  four  different  methods, 
including  combinations  from  two  different  Pn-wave  windows  (a  0.7  km/s  group-velocity 
window  around  the  IASP91  arrival  times  and  a  fixed  4-second  window  after  the  first 
arriving  P  wave)  and  two  different  displacement  components  (vertical  component  and  the 
displacement  rotated  to  the  Pn  incident  direction).  The  methods  generally  exhibit 
consistent  results.  Therefore,  we  only  present  the  result  from  the  vertical  component  and  a 
0.7  km/s  group-velocity  window. 

7  Conclusions 

We  present  a  method  to  separate  the  geometric  spreading  and  attenuation  from  the 
seismic  Pn-wave  data.  The  frequency-dependent,  log-quadratic  spreading  function  of  Yang 
et  al.  [2007]  and  the  power-law  Q  model  are  adopted  for  Pn  propagation.  A  high-accuracy 
broadband  Pn-wave  dataset  from  the  recent  North  Korean  nuclear  explosions  is  used  to 
constrain  the  model  parameters.  The  geometric  spreading  function  and  regional  Pn  Q  are 
obtained  for  Northeast  China  and  the  Korean  Peninsula.  By  taking  the  two-station 
amplitude  ratios  and  correcting  for  the  known  spreading  function,  we  remove  the  effects 
of  the  source  and  crust  legs  from  Pn  data,  obtaining  the  P-wave  attenuation  information 
along  the  pure  upper  mantle  path.  Combining  both  the  single-  and  two-station  data,  the 
upper  mantle  P-wave  attenuation  distribution  is  obtained  using  a  formal  tomographic 
approach.  The  current  method  can  be  applied  to  earthquake  data  as  well  as  Sn  waves.  We 
compared  the  Pn  waves  across  China  and  through  the  Japan  Sea.  The  results  reveal 
prominent  differences  for  Pn  waves  crossing  different  paths,  with  the  high-frequency  Pn 
signal  propagating  more  efficiently  through  the  continental  path  compared  with  the 
oceanic  path. 
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Table  1*  Event  Parameters  Used  in  This  Study 


North  Korean 
Nuclear  Test 

Date 

(yyyy/mm/dd) 

Time 

(UTQ 

Latitude 

(*N) 

Longitude 

TE) 

Depth 

(km) 

mb 

(USGS) 

NKT1 

2006/10/09 

0135:28.00 

41.287 

129.108 

0.0 

4.2 

NKT2 

2009/05/25 

00:54:43.1 1 

41.294 

129.077 

0.0 

4.7 

NKT3 

2013/02/13 

02:57:51.27 

41.292 

129.073 

0,0 

5.1 

Table  2,  Model  Parameters  of  Pn  Spectral  Amplitude  Used  in  This  Study 

Model  Space  and  Prior  Information 

Best  Fit  Inverted  Model 

Parameter 

Mame 

Prior  Information 

Data  Range 

Reference 

(2.0-10.0  Hz) 

Source 

Seismic  moment  for 

201 3/02/1 3 

2.0E  + 16 

5.0E  +  15<Mo<1,GE  +  17 

Taylor  et  at,  [2002] 

(7.1 33  ±  1.409)E  +■  15 

2009/05/25 

3.0E+  15 

6.GE  +  14<Mo<5.0E  +  15 

(4.015±0.591)E+15 

2006/10/09 

4.QE+14 

1.0E  +  14<Mo<6.0E  +  14 

(5.080  ±0,663)E  +  14 

b 

Slope  for  logMo'^g/c 

4.73 

3,5<b<5,5 

Xie  and  Patton  [1999] 

5,1 12  ±0329 

relation 

Geometric  Spreading 

4u 

Coefficient  for 

Pn  spreading 

-0217 

-0,239  <  In  <-0.195 

Yang  et  ol.  [2007] 

-0.213  ±0,010 

4n 

1,79 

1.611  <f12<  1969 

1.789  +  0.089 

£13 

3,16 

2,844  <£u  <3.476 

3.1 14  ±0,131 

-1.94 

-2,134<£21  <—1.746 

-1.993  ±0,102 

422 

8.43 

7,587  <^2^5.273 

8.459+0.378 

18.6 

1 6.74  <f2J<  20.46 

18.328+ 0.625 

4%  i 

-3.39 

-3,279<^  <-3051 

-3.407  ±0.183 

^32 

9.94 

8-946 <^32  <10,934 

9,787  ±0.509 

^33 

20,7 

18.63  <  £3  j  <22.77 

20.527  ±0.91 6 

Apparent  Q  Model 

Qq 

1  Hz  Q 

204-325 

200  <  Otj  <  600 

e.g,r  Xie  [2007]  and 

237  (194  -  290) 

J? 

Frequency 

dependence 

0.14-0.74 

— 0l2<j/<1.2 

Yang  [201 1] 

0.361  +0.071 

100 
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Figure  1.  Map  depicting  locations  of  the  North  Korean  test  site  (red  star  labeled  NKTS) 
and  seismic  stations  used  in  the  study,  including  CNDSN  (upset  down  triangles),  GSN 
(solid  circles),  SASK  (squares),  NECsaids  (open  circles),  and  F-NET  (triangles).  The  large 
circle  indicates  the  1000-km  epicentral  distance  from  the  NKTS,  with  the  pink  and  blue 
colors  marking  stations  with  continental  and  oceanic  paths,  respectively. 
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Figure  2.  (a)  A  sample  seismogram  recorded  at  station  LN.LYA,  (b)  its  normalized 
energy,  and  (c)  the  record  section  of  normalized  waveform  energy  for  all  data  used  in  this 
study.  Color  bar  indicates  the  normalized  energy  level.  The  regional  phases  Pn,  Pg,  and  Lg  are 
labeled  in  (b)  and  (c). 
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Figure  3.  Pn  wave  spectral  amplitude  measurement,  (a)  Vertical-component  Pn 
seismogram,  with  the  shaded  areas  indicating  windows  for  Pn  and  the  pre-signal  noise,  (b) 
Pn  phase  sampled  by  a  0.7-km/s  group  velocity  window,  (c)  Pn-wave  and  noise  spectra, 
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(d)  signal-to-noise  ratio  (SNR),  and  (e)  Pn  spectra  after  noise  correction.  Note  that  the  data 
samples  for  SNR  <  2.0  are  eliminated. 
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Figure  4.  Pn  spectral  amplitudes  versus  distance  at  selected  frequencies.  The  source 
excitation  spectra  have  been  removed.  Different  symbols  indicate  data  from  NKT1 
(triangles),  NKT2  (circles),  and  NKT3  (crosses).  The  solid  blue  lines  are  the  best-fit  Pn 
spectra  with  data  between  2.0  and  10.0  Hz,  and  the  shaded  areas  represent  their  standard 
deviations.  The  frequency  (f),  number  of  samples  (N),  apparent  Pn  Q,  and  correlation 
coefficient  (R)  are  labeled  in  each  panel. 
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Figure  5.  The  inverted  Pn  models,  (a)  and  (b)  The  Pn  geometric  spreading  functions 
versus  distance  at  2. 1  and  9. 1  Hz.  (c)  and  (d)  The  Pn  geometric  spreading  functions  versus 
frequency  at  500  and  1,200  km.  (e)  Source  excitation  functions  for  three  nuclear 
explosions,  and  (f)  regional  average  Pn  Q  versus  frequency.  The  solid  lines  are  the  result 
for  Northeast  China  and  the  Korean  Peninsula  obtained  using  2.0  to  10.0  Hz  data,  and  the 
dashed  lines  represent  the  model  given  by  Yang  et  al.  [2007]. 
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Figure  6.  Sketch  depicting  the  Pn-wave  propagation  path  from  the  NKTS  to  LN.LYA. 
The  crust  model  is  based  on  CRUST  1.0. 
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Figure  7.  (a)  -  (c)  Spreading-corrected  Pn  spectral  ratios  (light  gray  crosses)  versus 

interstation  distance  at  3.0,  5.2  and  9.1  Hz.  The  solid  lines  indicate  linear  regressions,  and 

their  slopes  represent  the  regional  average  uppermost  mantle  P-wave  Q.  (d )  Regional 

average  uppermost  mantle  P-wave  Q  (circles)  versus  frequencies.  The  solid  line  is  the 

best-fit  power-law  Q  model  using  data  between  2.0  and  10.0  Hz. 
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Figure  8.  Results  for  the  upper  mantle  P-wave  Q  model,  (a)  Inverted  Q  map  at  3.0  Hz,  (b) 
the  broadband  Q  obtained  by  averaging  logarithmic  Q  between  2.0  and  10.0  Hz,  (c)  3  Hz 
ray-path  coverage,  in  which  the  blue  and  pink  lines  indicate  the  single-  and  two-station 
paths,  respectively,  and  (d)  2.0°  x  2.0°  checkerboard  test  for  3  Hz  resolution  check. 
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Figure  9.  Pn  spectral  amplitudes  at  distance  1,000  km  versus  the  azimuth.  The  triangles, 
circles  and  crosses  represent  NKT1,  NKT2  and  NKT3,  respectively.  Black,  blue  and  red 
colors  indicate  0.8-,  7-  and  10.0-Hz  data,  respectively.  Solid  circles  with  error  bars 
represent  their  mean  values  and  standard  deviations,  obtained  within  30-degree  azimuth 
windows.  The  Pn  source  excitation  functions  are  removed  from  the  data. 
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